\section{Introduction}

A famous conjecture of  Hadwiger~\cite{H43} states that every graph $G$ of chromatic number $\chi (G)\ge t$ contains  $K_{t}$ as a minor. The conjecture is only known to be true for $t\le 6$, and  probably hard for larger values of $t$, as the proofs for cases $t=5, 6$ already depend on the Four Color Theorem. In general, results of Kostochka~\cite{K84} and, independently, Thomason~\cite{T84} imply that every graph~$G$ with $\chi(G)\ge t$ contains  a clique on $\lceil \frac {t}{f(t)} \rceil$ vertices as a minor,  where $f(t)\in O(\sqrt{\log t})$. Recent breakthroughs of Norin, Postle, and Song~\cite{NPZ20} and of Postle~\cite{P20} have improved on this general bound, but it is still unknown whether $f(t)=c$ for some constant $c$.

It is easy to see that $\chi(G)\ge \lceil\frac n{\alpha(G)}\rceil$ for every $n$-vertex graph $G$, where $\alpha(G)$ denotes the independence number of~$G$. So,
Hadwiger's conjecture, if true, implies that every $n$-vertex graph contains a clique minor of order at least $\lceil\frac n{\alpha(G)}\rceil$. Duchet and Meyniel~\cite{DM82} conjectured that this holds. Note that the order of this clique minor would be
 best possible as $G$ could be the disjoint union of cliques.

Providing evidence for their conjecture,
Duchet and Mey\-niel~\cite{DM82} proved that every graph~$G$ on $n$ vertices contains a clique minor of order $\lceil\frac{n}{2\alpha(G)-1}\rceil$. There have been several improvements on the order of the clique minor~\cite{BLW11,F10,KPT05,KS07,W07}, the best bound, due to Balogh and Kostochka~\cite{BK11}, being $\lceil\frac{n}{c\alpha (G)}\rceil$, where $c$ is a constant with $c<1.95$.

The focus of this paper is an analogous result replacing minors  with  immersions.
A graph~$G$ is said to contain another graph $H$ as an \emph{immersion} if there exists an injective function $\phi\colon V(H)\rightarrow V(G)$ such that:
\begin{enumerate}[(I)]
\item For every $uv\in E(H)$, there is a path in $G$, denoted $P_{uv}$, with endpoints $\phi(u)$ and~$\phi(v)$.
\item The paths in $\{P_{uv} \mid uv\in E(H) \}$ are pairwise edge disjoint.
\end{enumerate}
The vertices in $\phi(V(H))$ are called the \emph{branch vertices} of the immersion. If branch vertices are not allowed to appear as interior vertices on paths $P_{uv}$, the immersion is called {\em strong}.

The minor order and the immersion order are not comparable. The class of planar graphs, while excluding~$K_5$ as a minor, contains all cliques as immersions. On the other hand, the class of graphs with maximum degree at most $d$, while excluding $K_{d+2}$ as an immersion, contains all cliques as minors.
However, the two relations do share some important similarities. Both of them are well-quasi-orders~\cite{RS04,RS10}, and both notions are weakenings of the topological minor relation. Possibly inspired by such similarities, Lescure and Meyniel~\cite{LM89}  proposed an analogue of Hadwiger's conjecture for strong immersions. Later, Abu-Khzam and Langston~\cite{AKL03} weakened their conjecture to the following form.

\begin{conjecture}[\cite{AKL03, LM89}] \label{hadimm}
Every graph $G$ with $\chi(G)\ge t$ contains $K_{t}$ as an immersion.
\end{conjecture}

This conjecture and its strong version have received much attention recently, and have been tackled with more success than their minor counterpart. The cases $1\le t \le 4$  follow from the fact that Haj\'os' Conjecture
 is true for these cases~\cite{D52}. The cases $5\le t \le 7$ of Conjecture~\ref{hadimm} were established by Lescure and Meyniel~\cite{LM89} and by DeVos, Kawarabayashi, Mohar, and Okamura~\cite{Detal10}.

For general values of $t$, the first linear lower bound for Conjecture~\ref{hadimm} was given by DeVos, Dvo\v{r}\'ak, Fox, McDonald, Mohar, and Scheide~\cite{Detal14}. They proved that every graph~$G$ contains an immersion of a clique on $\bigl\lceil \frac{\chi (G)}{200} \bigr\rceil$ vertices. Dvo\v{r}\'ak and Yepremyan~\cite{DY17} improved this bound to $\bigl\lceil \frac{\chi (G)-7}{11} \bigr\rceil$. The best currently known lower bound for Conjecture~\ref{hadimm} is due to Gauthier, Le and Wollan~\cite{GLW17} who showed that every graph $G$ contains an immersion of a clique on $\bigl\lceil \frac{\chi (G)-4}{3.54} \bigr\rceil$.  This implies the following.

\begin{theorem}[Gauthier, Le and Wollan~\cite{GLW17}]\label{wow}
Every $n$-vertex graph $G$ contains  an immersion of a clique on  $\bigl\lceil \frac{n}{3.54\alpha(G)}-1.13\bigr\rceil$ vertices.
\end{theorem}

It seems natural to ask whether the bound from Theorem~\ref{wow} can be improved without necessarily improving Gauthier, Le and Wollan's underlying result that relates the chromatic number and the size of clique immersions,  as Duchet and Meyniel did in the context of graph minors.
The first attempt in this direction (actually earlier than~\cite{GLW17}) has been carried out by Vergara~\cite{V17}. She conjectured that
 every $n$-vertex graph with  independence number 2 contains an immersion of $K_{\lceil \frac n2\rceil}$, and showed that this conjecture is equivalent to Conjecture~\ref{hadimm} for graphs of independence number 2.
 In support of her conjecture, Vergara  proved that every graph on $n$ vertices and independence number 2  contains a $K_{\lceil \frac n3\rceil}$-immersion. Her result was improved by Gauthier, Le and Wollan~\cite{GLW17} as follows.

\begin{theorem}[Gauthier, Le and Wollan~\cite{GLW17}]\label{dos}
Every $n$-vertex  graph with independence number~$2$ contains $K_{2\lfloor \frac n5 \rfloor}$ as an immersion.
\end{theorem}

Additionally, Vergara's conjecture has been proved to hold for graphs with small forbidden subgraphs~\cite{Q20}.
Extending the conjecture of Vergara, we propose the following conjecture for graphs of arbitrary independence number.

\begin{conjecture}\label{ourconj}
Every $n$-vertex  graph $G$ contains $K_{\lceil \frac n{\alpha(G)}\rceil}$ as an immersion.
\end{conjecture}

This conjecture is best possible, since (as mentioned earlier for the case of minors) the graph  $G$ could be the disjoint union of cliques.

Our main result improves on Theorem~\ref{wow} for every $G$ with $\alpha(G)\ge 3$, giving more evidence in support of Conjecture~\ref{ourconj}, and thus of Conjecture~\ref{hadimm}.

\begin{theorem}\label{all}
Let $G$ a graph on $n$ vertices  with independence number  $\alpha\in\mathbb N$. Then $G$ contains an immersion of a clique on at least $\lfloor  \frac {n}{2.25\alpha  -f(\alpha)} \rfloor - 1$ vertices, where $f(3)=f(4)=2.25$ and $f(\alpha)=2.25+\sum_{i=6}^{\alpha+1} \frac 1{2i}$, if $\alpha\ge 5$.
\end{theorem}

In particular, Theorem~\ref{all}  implies that any $n$-vertex graph with independence number $\alpha\le 10$ contains an immersion of a clique on $\lfloor \frac {n}{2 \alpha} \rfloor - 1$ vertices.

In the literature, several other types of immersions have been distinguished (apart from immersions and strong immersions). If all the paths~$P_{uv}$ from~(I) have length at most $k$, then~$H$ is called a \emph{$k$-immersion} of $G$. If  all the paths $P_{uv}$ have odd length then~$H$ is an \emph{odd immersion}.

Actually, Theorems~\ref{wow} and~\ref{dos} as well as   the precursor results give  strong immersions. Moreover, one can read from
the proof of Theorem~\ref{dos}  that the immersion from that theorem is a strong odd $5$-immersion. While no such information can be deduced for Theorem~\ref{wow}, it is not difficult to verify that the proof of our Theorem~\ref{all} does imply a stronger statement.

\begin{remark}\label{more}
The immersion from Theorem~\ref{all} is  a strong odd $(2\alpha -1)$-immersion.
\end{remark}

In particular, this means that our result not only gives further evidence for Conjectures~\ref{hadimm} and~\ref{ourconj}, but also for a stronger conjecture of Churchley~\cite{C17} stating that every graph of chromatic number at least~$t$  contains $K_{t}$ as an odd immersion.

The structure of the paper is as follows. In the next section (Section~\ref{secaalpha}) we introduce the auxiliary graphs $A_\alpha$ and prove a crucial lemma, Lemma~\ref{parti}. This lemma will be used at the end of the proof of Theorem~\ref{all}, which will be given  in Section~\ref{secall}.

An extended abstract~\cite{nos} announcing a weaker version of Theorem~\ref{all}  includes an alternative proof for the case $\alpha(G)=3$.

\section{The graph $A_\alpha$ and its minimal cuts}\label{secaalpha}

In this section we give key definitions and prove some crucial ingredients for the proof of Theorem~\ref{all}.

\subsection{Minimal cuts and blow-ups}

We first need to introduce some standard notation on (vertex-)cuts. Let $G$ be a graph. For $A,B\subseteq V(G)$, an \emph{$A$-$B$ cut} is a set $C\subseteq V(G)$ such that $G\setminus C$ contains no path from $A\setminus C$ to $B\setminus C$. (Slightly abusing notation, we might sometimes write $a$-$B$ cut if $A=\{a\}$.)
If $C$ is an $A$-$B$ cut we also say $C$ is \emph{a cut separating $A$ and $B$}. We say $C$ is \emph{minimal} if for every $x\in C$, we have that $C\setminus\{ x\}$ is not an $A$-$B$ cut. We say $C$ is \emph{minimum} if $|C|\le |C'|$ for every other $A$-$B$ cut $C'$.

Let $\mathbb{N}$ denote the set of non-negative integers. For a graph $G$ and $f\colon V(G)\rightarrow\mathbb{N}$ a function, an \emph{$f$-blow-up} of $G$ is a graph that can be obtained from $G$ by replacing each vertex $v\in V(G)$ with an independent set $B(v)$ satisfying $|B(v)|=f(v)$, and $B(v)\cap B(u)=\varnothing$ if $v\ne u$,  and having the edge $xy$, for $x\in B(v), y\in B(u)$, if and only if $uv\in E(G)$.  We  let $B(X)=\bigcup_{v\in X}B(v)$ for every $X\subseteq V(G)$.

We leave the proof of the following lemma to the reader.
\begin{lemma}\label{blowcut}
Let $G$ be a graph, let $f\colon V(G)\rightarrow \mathbb{N}$, and let $X, Y\subseteq V(G)$. Let $B$ be an $f$-blow-up of~$G$, and let $C_B$ a minimum $B(Y)$-$B(X)$ cut in $B$.  Then $C:=\{v\in V(G) : B(v)\cap C_B\neq\varnothing\}$ is a minimal $X$-$Y$ cut in $G$  and $C_B=\bigcup_{v\in C}B(v)$.
\end{lemma}

\subsection{The  auxiliary graph $A_\alpha$}

We now define a family of graphs $\{ A_\alpha \}_{\alpha\ge 2}$. These graphs will represent, in a simplified way, the structure that arises when trying to find the desired immersion for our theorem by induction.

There is another auxiliary graph involved in our proof, and this is a blow-up of~$A_\alpha$. The function which defines this blow-up varies with the graph $G$, this blow-up being the connection between our arbitrary graph $G$ (of independence number $\alpha$) and our fixed graph $A_\alpha$. So it is good to think of each vertex (except $d_\alpha$) of the auxiliary graph $A_\alpha$ as representing a subsets of $V(G)$, and each edge representing a complete bipartite subgraph. Yet, the blow-up of $A_\alpha$ is not necessarily (isomorphic to) a subgraph of $G$. The exact relation between the blow-up of $A_\alpha$ and $G$ will be clarified later. For the final structural analysis however, it is easier (and sufficient) to restrict our attention to~$A_\alpha$.

For a given $\alpha\ge 2$ we define the graph $A_\alpha$ as the one having the following as its vertices,
\begin{itemize}\itemsep -2pt
\item a vertex $d_\alpha$;

\item  a vertex $d_S$ for every nonempty $S\subseteq \{1,\dots ,\alpha-1\}$; and

\item a vertex $x_S$ for every $S \subseteq \{1,\dots ,\alpha\}$ with $|S|\ge 2$;
\end{itemize}

\noindent
and the following as its edges
\begin{itemize}\itemsep -2pt
\item the edge $d_\alpha x_S$, for every $S \subseteq \{1,\dots ,\alpha\}$ with $\alpha\in S$,  and $|S|\ge 2$;

\item the edge $d_Sx_T$, for each pair $S,T\subseteq \{1,\dots ,\alpha\}$ with $\alpha\notin S$, $|T|\ge 2$, and $S\cap T\ne\varnothing$;

\item the edge $x_Sx_T$, for each pair $S,T\subseteq \{1,\dots ,\alpha\}$ with $|S|, |T|\ge 2$, and $S\cap T\ne\varnothing$.
\end{itemize}

We write $$D_\alpha:=\{d_S \mid S \subseteq \{1,\dots ,\alpha-1\}, S\ne \varnothing\},$$ and observe that
\begin{equation}\label{neighbs}
N_{A_\alpha}(d_T)\subseteq N_{A_\alpha}(d_S)
\end{equation}
for all $T\subset S$, with $d_T, d_S \in D_\alpha$.

Since we plan to use Menger's theorem later for finding the connections between the new branch vertex and the old ones, we will be interested in minimum cuts of the original graph~$G$. Here we investigate the minimal cuts of the graphs $A_\alpha$, and later use Lemma~\ref{blowcut} to translate what we find here into information about the cuts of $G$. The following lemma contains all the important properties of minimal cuts we will use.

\begin{lemma}\label{Aalpha}
Let $\alpha\ge 2$ and let $C$ be a minimal $d_\alpha$-$D_\alpha$ cut in $A_\alpha$ such that $C\ne \{d_\alpha\}$. Then all of the following hold, for each $S\subset \{1,\dots ,\alpha-1\}$.
\begin{enumerate}[(a)]
\item\label{aalphad}
If  $d_S\notin C$, then $d_T\notin C$ for every nonempty $T\subset S$.
\item\label{aalphax}
If  $x_{S\cup\{\alpha\}}\notin C$, then $x_{T\cup\{\alpha\}}\notin C$ for every nonempty $T\subset S$.
\item\label{aalpha1}
If $d_S\notin C$, and if $x_{T\cup\{\alpha\}}\notin C$ for some nonempty $T\subseteq \{1,\dots ,\alpha-1\}$, then $x_R \in C$ for every $R\subseteq \{1,\dots ,\alpha-1\}$ with $R\cap S\ne\varnothing\ne R\cap T$ and $|R|\ge 2$.
\item
\label{aalpha2}
We have $d_S\notin C$ if and only if $x_{T\cup\{\alpha\}}\in C$ for every nonempty $T\subseteq \{1,\dots ,\alpha-1\}$ with $T\cap S\ne \varnothing$.
\end{enumerate}
\end{lemma}

\begin{proof}
We first prove~\eqref{aalphad}. Assume there is $T\subset S$ with $d_T\in C$. Since $C$ is minimal, $A_\alpha\setminus (C-d_T)$ contains a $d_T$-$d_\alpha$ path. But then, by~\eqref{neighbs}, either  $A_\alpha\setminus C$ contains a $d_S$-$d_\alpha$ path, which is impossible, or $d_S\in C$, which is as desired. This proves~\eqref{aalphad}.

Property~\eqref{aalphax} can be proved using similar arguments as for~\eqref{aalphad} since $N_{A_\alpha}(x_{T\cup\{\alpha\}})\subseteq N_{A_\alpha}(x_{S\cup\{\alpha\}})$. And property~\eqref{aalpha1} holds since $d_Sx_Rx_{T\cup \{\alpha\}}d_\alpha$ is a path.

Let us now prove~\eqref{aalpha2}.
Suppose first that we have $d_S\notin C$. If $T\cap S\ne \varnothing$, consider the path $d_Sx_{T\cup\{\alpha\}}d_\alpha$. This path has to meet $C$, and thus we have  $x_{T\cup\{\alpha\}}\in C$.

Suppose now that $x_{T\cup\{\alpha\}}\in C$ for every nonempty $T$ with $T\cap S\ne \varnothing$. For a contradiction, assume $d_S\in C$. Since~$C$ is minimal, there must be a $d_S$-$d_\alpha$ path in $A_\alpha\setminus(C-d_S)$. Let $P=d_Sx_{R_k}\dots x_{R_1}x_{R_0}d_\alpha$ be a minimum length path among these. By hypothesis we have that $k\ge 1$.

Let us show that there exists $U\neq\varnothing$ such that $U\subseteq R_k$ and $x_{U\cup\{\alpha\}}\in C$. If $k=1$, then by definition of $P$, we have $S\cap R_1\ne \varnothing$. So we can take $U=S\cap R_1$ (notice that $x_{U\cup \{\alpha\}}\in C$ by hypothesis). If $k\ge 2$, we note that since $P$ is of minimum length, the path $d_Sx_{R_k}x_{R_k\cup\{\alpha\}}d_\alpha$ meets $C$, which implies $x_{R_k\cup\{\alpha\}}\in C$. So we can take $U=R_k$.

By minimality of $C$, the graph $A_\alpha\setminus(C-x_{U\cup\{\alpha\}})$ contains a $D_\alpha$-$d_\alpha$ path that passes through $x_{U\cup\{\alpha\}}$. Let $x_W$, or alternatively $d_W$, be the vertex before $x_{U\cup\{\alpha\}}$ on this path, and $P'$ the subpath of this path that joins $D_\alpha$ with either $x_W$ or $d_W$. If we had $\alpha\in W$, the vertex before $x_{U\cup\{\alpha\}}$ would not be $d_W$, but rather $x_W$ and adjacent to $d_\alpha$, and there would be a $D_\alpha$-$d_\alpha$ path that avoids $C$, a contradiction. Then we have $\alpha\notin W$ and, since $U\subseteq R_k$, we have  $\varnothing \ne W\cap U\subseteq W\cap R_k$. Then  $P'x_{R_k}\dots x_{R_1}x_{R_0}d_\alpha$ is or contains a $D_\alpha$-$d_\alpha$ path that avoids $C$, a contradiction.
\end{proof}

We are now ready to show the crucial lemma for the proof of Theorem~\ref{all}.

\begin{lemma}\label{parti}
Let $\alpha\ge 2$ and let $C$ be a minimal $d_\alpha$-$D_\alpha$ cut in $A_\alpha$ such that $C\ne \{d_\alpha\}$. Then there are two (possibly empty) disjoint sets $J_1$ and $J_2$, such that $J_1\cup J_2=\{1,\dots ,\alpha -1\}$ and
\begin{enumerate}[(i)]
\item $x_S\in C$ for every $S\subseteq \{1,\dots ,\alpha -1\}$ with $J_1\cap S\neq\varnothing\neq J_2\cap S$;
\item $x_{S\cup\{\alpha\}}\in C$ for every $S\subseteq \{1,\dots ,\alpha -1\}$ with $J_2\cap S\neq\varnothing$;
\item\label{vacio} If $J_2=\varnothing$ then $d_S\in C$ for every nonempty $S\subseteq \{1,\dots ,\alpha -1\}$.
\end{enumerate}
\end{lemma}

\begin{proof}
Let $J_1$ be the smallest subset of $\{1,\dots ,\alpha -1\}$ such that  for every nonempty $T$ with $x_{T\cup \{\alpha\}}\notin C$ we have $T\subseteq J_1$. Let $J_2$ be the smallest subset of $\{1,\dots ,\alpha -1\}$ such that  for every $S\subseteq \{1,\dots , \alpha -1\}$ with $d_S \notin C$ we have $S\subseteq J_2$. Clearly \emph{(iii)} holds for this choice.

By Lemma~\ref{Aalpha}~\eqref{aalpha2}, we know that if $d_S \notin C$ then $x_{T\cup \{\alpha\}}\in C$ for every $S,T\subseteq \{1,\dots ,\alpha -1\}$ with $S\cap T\ne \varnothing$. By the choice of $J_1$ and $J_2$ as minimal, this implies that $J_1\cap J_2=\varnothing$.
Moreover, if there exists $i\in \{1,\dots ,\alpha -1\}\setminus (J_1\cup J_2)$, then we have  $x_{\{ i\}\cup\{\alpha\}}\in C$ and $d_{\{ i\}}\in C$, which contradicts Lemma~\ref{Aalpha}~\eqref{aalpha2}. Thus, we also have  $J_1\cup J_2=\{1,\dots ,\alpha -1\}$.

By Lemma~\ref{Aalpha}~\eqref{aalphad} and~\eqref{aalphax}, it follows that for every $i\in J_1$ and every $j\in J_2$ we have $x_{\{i\}\cup\{\alpha\}}\notin C$ and $d_{\{ j\}}\notin C$. This implies, by Lemma~\ref{Aalpha}~\eqref{aalpha1}, that for every $i\in J_1$ and every $j\in J_2$ we have $x_S\in C$ for every $S\subseteq \{1,\dots ,\alpha -1\}$ with $i,j\in S$. Thus \emph{(i)} holds. Moreover,   Lemma~\ref{Aalpha}~\eqref{aalpha2} implies that for every $j\in J_2$ we have $x_{S\cup\{\alpha\}}\in C$ for every $S\subseteq \{1,\dots ,\alpha -1\}$ with $j\in S$, and therefore, also \emph{(ii)} holds.
\end{proof}

\section{The proof of Theorem~\ref{all}}\label{secall}

We first need the following elementary lemma.

\begin{lemma}\label{elem}
Let $\beta_1=1$, $\beta_2=\frac{5}2$, $\beta_3=\frac 92$, $\beta_4=\frac {27}4$ and $\beta_i=\frac94 (i-1) - \sum_{j=6} ^{i+1} \frac 1{2j}$  for $i\ge 5$. Let~$s,t$ and $i\ge 4$ be integers with $s+t=i$, and $2\le s\le t\le i-2$. Then we have $$\beta_i\ge  \beta_s+\beta_t +\frac{2i -1}i.$$
\end{lemma}
\begin{proof}
If $i=4$, then $s=t=2$ and the result holds with equality. So we assume $i\ge 5$.  If $s=2$, notice that we have $$\beta_{i}=\beta_{i-2}+\frac92 -\frac 1{2i}- \frac 1{2(i +1) } > \beta_2 +\beta_{i-2} +\frac{2i -1}i ,$$ as desired. When $s=3$ we must have $i\ge 6$, and thus $$\beta_{i}\ge \beta_{i-3}+\frac{27}4 - \frac 1{10} -\frac 1{2i } - \frac 1{2(i +1) } > \beta_3 +\beta_{i-3} +\frac{2i -1}i .$$The case $s=4$ follows analogously. When $s\ge 5$, it is not hard to see that we have
\begin{equation*}
\beta_{i} = \beta_{t}+\frac {9s}4- \sum_{j=t+2}^{i+1}\frac 1{2(i +1)} >  \beta_t+  \beta_s + \frac{2i -1}i ,
\end{equation*}
which is as desired.
\end{proof}

Theorem~\ref{all} follows immediately from Lemma~\ref{elem} together with the following result.
\
 \begin{theorem}\label{thmX}
 Let  $\beta_1,\beta_2,\beta_3, \ldots$ be rationals such that $\beta_1=1,\beta_2=\frac 52$, while for all $i\ge 3$ \\
\phantom{yap} $\bullet$ \, $\beta_i\ge \beta_{i-1}+2$ and \\
\phantom{yap} $\bullet$ \,  for all  $s, t\in\{2, \ldots, i-2\}$ with $s+t=i$  we have $\beta_i\ge \beta_s+\beta_t +\frac{2i-1}i.$ \\
Let $G$ be a graph on~$n$ vertices with independence number at most $\alpha$. Then $G$ contains an immersion of a clique on at least $\lfloor \frac n{\beta_\alpha} \rfloor - 1$ vertices.
\end{theorem}

\begin{proof}[Proof of Theorem~\ref{thmX}]
We proceed by  induction on $\alpha$. The statement is trivially true for $\alpha=1$, and true for $\alpha=2$ by Theorem~\ref{dos}. We now prove it for~$\alpha\ge 3$. In turn, we do this by induction on $n$. It is easy to see that the statement holds for all $1\le n<3\beta_\alpha$ as we then only need to find a $K_1$-immersion.

Note that we can assume that $\alpha (G)=\alpha$, and consider a maximum independent set $I=\{a_1, \dots, a_\alpha\}$ of~$G$. By  induction on $n$, we know that $G-I$ contains an immersion  of a clique on $\lfloor \frac{n-\alpha}{\beta_{\alpha(G-I)}} \rfloor - 1\ge \lfloor \frac{n-\alpha}{\beta_{\alpha}} \rfloor - 1$ vertices, with set of branch vertices~$M$. Our plan is to  add a new branch vertex to this immersion, or, failing this, to find  some other immersion of  a clique on at least $\lfloor \frac n{\beta_\alpha} \rfloor - 1$ vertices in $G$. If we succeed in adding a new branch vertex to the immersion, we will have constructed an immersion on $\lfloor \frac{n-\alpha}{\beta_{\alpha}} \rfloor - 1+1\ge \lfloor \frac{n}{\beta_{\alpha}} \rfloor - 1$ vertices, which is as desired.

Set $Q:=V(G)\setminus (I\cup M)$ and note that
\begin{equation}\label{Qall}
|Q|\ = \ n-\alpha -(\Big\lfloor \frac{n-\alpha}{\beta_\alpha} \Big\rfloor - 1)\ > \ \frac{\beta_\alpha-1}{\beta_\alpha}n - \alpha.
\end{equation}

Set $N_i:=N(a_i)$ and $\bar N_i:= V(G)\setminus (N(a_i)\cup \{a_i\} )$ for $i=1,\dots ,\alpha$. Since $G$ has independence number~$\alpha$, we know that if $S\subseteq\{ 1,\dots,\alpha\}$ then $G[\bigcap_{i\in S}\bar N_i]$ has independence number at most $\alpha-|S|$. Hence, we may assume that
\begin{equation}\label{nonneighS}
\Bigl| \bigcap_{i\in S} \bar N_i\Bigr| <\frac{\beta_{\alpha-|S|}}{\beta_\alpha}\, n, \mbox{ for every } S\subseteq\{ 1,\dots,\alpha\}, \, 1\le|S|\le \alpha -1
\end{equation}
as otherwise, by induction on $\alpha$, there would be some $S$ such that $G[\bigcap_{i\in S} \bar N_i]$ contains an immersion of $K_{\lfloor \frac{n}{\beta_{\alpha}} \rfloor -1}$.

Let us prove next that we have
\begin{equation}\label{neeew}
\text{$|  M\cap \bar N_i | < | \bigcup_{j \neq i} Q\cap N_i \cap N_j|$ for each $i \in \{1,\dots ,\alpha \}$.}
\end{equation}
Indeed for a contradiction, suppose there exists $i$ such that $| M\cap \bar N_i  | \geq | \bigcup_{j \neq i} Q\cap N_i \cap N_j |$. Notice that $Q = (\bigcap_{j \neq i} Q\cap \bar N_j) \cup ( \bigcup_{j\neq i} Q\cap N_i \cap N_j ) \cup ( Q\cap \bar N_i )$. Since $\bigcap_{j \neq i}\bar  N_j$ is a clique, we obtain
\begin{equation*}
\begin{split}
|Q| &\leq \Big\lfloor \frac{n}{\beta_\alpha} \Big\rfloor -2 +| M\cap \bar N_i  |+ |Q\cap \bar N_i |\\
 &\leq \frac{n}{\beta_\alpha}  -2 + |\bar N_i\setminus I| \\
&< \frac{n}{\beta_\alpha} -\alpha -1 + \frac{\beta_{\alpha -1}}{\beta_\alpha}n \\
&= \frac{\beta_{\alpha-1}+1}{\beta_\alpha}n -\alpha -1 ,
\end{split}
\end{equation*}
where the third inequality comes from~\eqref{nonneighS}. This contradicts \eqref{Qall}, since $\alpha\ge 3$ and therefore by hypothesis we have $\beta_\alpha\ge \beta_{\alpha-1}+2$.

Without loss of generality we can assume $| M\cap \bar N_\alpha | =\min_{ i\in \{ 1,\dots ,\alpha \} } |M\cap \bar N_i|$ and so, by~\eqref{neeew},
\begin{equation}\label{boundapares}
| M\cap \bar N_\alpha  | < \Bigl|  \bigcup_{j \neq i} Q\cap N_i \cap N_j\Bigr|, \mbox{ for every } 1\le i\le \alpha.
\end{equation}

This assumption also implies
\begin{equation}\label{pendiente}
|M\cap N_\alpha | \ge \frac{1}{\alpha}|M|,
\end{equation}
since every vertex of $M$ is a neighbor of some vertex in $I$.

 We  claim that if
 $I_1$ and $I_2\ne \varnothing$ are disjoint and satisfy $I_1\cup I_2=\{1,\dots, \alpha-1\}$, then
 \begin{equation}\label{claimst}
 \Bigl| Q\cap \Bigl(\bigcup_{i\in I_1\cup \{\alpha\}} N_i\Bigr)\cap \Bigl(\bigcup_{j\in I_2} N_j\Bigr)\Bigr|\ge |M\cap \bar N_\alpha  |.
 \end{equation}
For the case $I_1=\varnothing$ and for the case $|I_2|=1$, \eqref{claimst} follows from \eqref{boundapares}. So we assume $I_1\ne \varnothing$ and $|I_2|\ge2$, and note that we have
\begin{equation*}
\begin{split}
\Bigl| Q\cap \Bigl(\bigcup_{i\in I_1\cup \{\alpha\}} N_i\Bigr)\cap \Bigl(\bigcup_{j\in I_2} N_j\Bigr) \Bigr| &\ge n - |M|- \Bigl|\bigcap_{i\in I_1\cup \{\alpha\}}\bar N_i\Bigr|- \Bigl|\bigcap_{j\in I_2}\bar N_j\Bigr|  \\
&> \Bigl(\beta_\alpha -1 - \beta_{\alpha -(|I_1|+1)} -\beta_{\alpha-|I_2|}\Bigr) \frac n{\beta_\alpha}\\
&\ge \frac{\alpha-1}\alpha\cdot\frac n{\beta_\alpha} \\ & > |M\cap \bar N_\alpha |,
\end{split}
\end{equation*}
where the second inequality comes from \eqref{nonneighS} and the fact that $|M|< \frac n{\beta_\alpha}$; the third inequality comes from the properties of the $\beta_i$; and the last inequality comes from~\eqref{pendiente} and (again) the fact that $|M|< \frac n{\beta_\alpha}$. Thus we have proved~\eqref{claimst}.

Let $\mathcal{P}_G$ be the set of all paths in $G$ starting in $a_\alpha$ and ending in $M\cap \bar N_\alpha $ which alternate between the set $M\cup Q$ and the set $I=\{a_1,\dots ,a_\alpha\}$ (i.e., each edge in such a path is incident to one vertex from each set) and have all internal vertices in $Q\cup I\setminus\{a_\alpha\}$\footnote{Note that if $P \in \mathcal{P}_G$ then the length of $P$ is odd and is at most $2\alpha-1$.}. Since these paths are all edge disjoint from the paths occupied by the immersion with branch vertices $M$, our goal is to show that there exists a collection of $|M \cap\bar N_\alpha |$ edge-disjoint paths in $\mathcal{P}_G$ such that all the paths in the collection have  different endpoints in $M\cap \bar N_\alpha $. Indeed, if we find such a collection, then the vertices in $M$ together with the vertex $a_\alpha$ would be  the branch vertices of an immersion of the desired size, and this would finish the proof.

Consider the function $f\colon V(A_\alpha)\rightarrow \mathbb{N}$ (where $A_\alpha$ is the graph defined in Section~\ref{secaalpha}) defined as follows\footnote{The value of $f(d_\alpha)$ does not really matter, as long as it is at least $|M\cap\bar N_\alpha|$.}:
\begin{itemize}\itemsep -2pt
\item $f(d_\alpha)= n$;

\item  $f(d_S)=|(M\cap \bigcap_{i\in S}N_i)\setminus \bigcup_{j\in \{1,\dots ,\alpha\}\setminus S}N_j|$, for every $S\subseteq \{1,\dots ,\alpha-1\}$;

\item $f(x_S)=|(Q\cap \bigcap_{i\in S}N_i)\setminus \bigcup_{j\in \{1,\dots ,\alpha -1\}\setminus S}N_j|$, for every $S \subseteq \{1,\dots ,\alpha\}$, with $|S|\ge 2$.
\end{itemize}

Let $B$ be an $f$-blow-up of $A_\alpha$. Let $\mathcal{P}_B$ be the set of all paths in $B$ with an endpoint in~$B(d_\alpha)$ and the other in $\bigcup_{S\subseteq\{1,\dots ,\alpha -1\}}B(d_S)$.

For $S,T$ with $S\ne\{\alpha\}$, let $M_S,Q_S,Q_T\subseteq V(G)\setminus I$ be the sets such that the definition of~$f$ gives $f(d_S)=|M_S|$, $f(x_S)=|Q_S|$ and $f(x_T)=|W_T|$. If $S\cap T\ne \varnothing$, then $Q_S\cup Q_T\subseteq N_i$ for every $i \in S\cap T$. With this in mind, we can see that how paths $\mathcal{P}_B$ represent paths in $\mathcal{P}_G$. For instance, a path going through the sets $B(d_\alpha),B(x_{\{ 2, \alpha\}}),B(x_{\{1,2\}}), B(d_{\{1\}})$ represents a path $a_\alpha q_{2\alpha} a_2 q_{12} a_1 m_1$ with $q_{2\alpha}\in Q_{\{2,\alpha\}}=(Q\cap N_2\cap N_\alpha)\setminus \bigcup_{j\in \{1,3,\dots ,\alpha -1\}}N_j$, $q_{12}\in Q_{\{1,\alpha\}}=(Q\cap N_1\cap N_2)\setminus \bigcup_{j\in \{3,4,\dots ,\alpha -1\}}N_j$ and $m_1\in M_{\{1\}} =(M\cap N_1)\setminus \bigcup_{j\in \{2,3,\dots ,\alpha\}}N_j$. Figure~\ref{theone}  illustrates  this example.

\begin{figure}[h]
 \centering
 \captionsetup{justification=centering}
 \bigskip
 \includegraphics[height=3 in]{immersion2.eps}
 \medskip
 \caption{Left: A graph $G$ with $\alpha(G)=3$, and a path $P\in \mathcal{P}_G$. Right: The path in $A_3$ corresponding to $P$ is shown with dashed edges.}
  \label{theone}
\end{figure}

Since the paths in $\mathcal{P}_G$ alternate between $M\cup Q$ and $I$, if $P, P'\in \mathcal{P}_G$ satisfy that $V(P)\setminus I$ and $V(P')\setminus I$ are disjoint, then we have that $P,P'$ are edge disjoint. Thus, if there is a family of $k$ vertex disjoint paths in $\mathcal{P}_B$, then  there is a family of $k$ edge-disjoint paths in $\mathcal{P}_G$ having different endpoints in $M\cap \bar N_\alpha $. So, if we can find a family of $|M\cap\bar N_\alpha|=\sum_{S\subseteq \{1,\dots, \alpha -1\}}f(d_S)$ vertex disjoint paths in $\mathcal{P}_B$, we are done.

By Menger's Theorem~\cite{M27}, it suffices to show that every (minimum) cut $C_B$ separating  $B(d_\alpha)$ from $\bigcup_{S\subseteq\{1,\dots , \alpha -1 \}}B(d_S)$ fulfills $|C_B|\ge |M\cap \bar N_\alpha|$. Lemma~\ref{blowcut}, in turn, tells us that $C_B=\bigcup_{v\in C}B(v)$, for some minimal cut $C$ separating $d_\alpha$ from $D_\alpha$ in $A_\alpha$. These facts and our choice of $f(d_\alpha)$ tell us that we can finish the proof by showing that

\begin{equation}\label{fini}
\text{$|\bigcup_{v\in C}B(v)|\ge |M\cap \bar N_\alpha|$ for every minimal $d_\alpha$-$D_\alpha$ cut $C$ with $C\ne \{d_\alpha\}$.}
\end{equation}

In order to see~\eqref{fini}, we take any such cut $C$ and consider the sets $J_1$ and $J_2$ given by Lemma~\ref{parti}. If $J_2\ne\varnothing$, then by definition of $f$ and later by~\eqref{claimst} we obtain
\begin{equation*}
\begin{split}
|\bigcup_{v\in C}B(v)| &\ge |\bigcup_{S\subseteq \{1,\dots ,\alpha -1\} : J_2\cap S\neq\varnothing}B(x_{S\cup\{\alpha\}})|+|\bigcup_{S\subseteq \{1,\dots ,\alpha -1\} : J_1\cap S\neq\varnothing\neq J_2\cap S}B(x_S)|\\
 &= \Bigl| Q \cap \Bigl(\bigcup_{i\in J_1\cup \{\alpha\}} N_i\Bigr)\cap \Bigl(\bigcup_{j\in J_2} N_j\Bigr)\Bigr| \\
&\ge |M \cap\bar N_\alpha|,
\end{split}
\end{equation*}
as desired. So we can assume $J_2=\varnothing$, in which case Lemma~\ref{parti}\emph{(\ref{vacio})} guarantees $d_S\in C$ for every $S\subseteq \{1,\dots ,\alpha -1\}$.  But then, the definition of $f$ gives $$|\bigcup_{v\in C}B(v)|\ \ge \ |\bigcup_{S\subseteq \{1,\dots ,\alpha -1\}}B(d_S)|\ =\ |M\cap\bar N_\alpha|.$$
This completes the proof of~\eqref{fini}, and thus the proof of the theorem.
\end{proof}

\section*{Acknowledgements}
We wish to thank two anonymous referees for their careful reading and suggestions.


\section{Introduction}\label{introduction}
The dynamical evolution of molecular clouds (MCs) in the interstellar medium (ISM) is determined by supersonic, compressible turbulent flows \citep{Larson1981,Solomon1987,Klessen2000,Heitsch2001,Ossenkopf2002,Elmegreen2004,Heyer2004,MacLow2004,Scalo2004,Krumholz2005,Ballesteros2007,McKee2007,RomanDuval2011,Padoan2014,Federrath2015}. The turbulent dynamics of the clouds plays a diverse and vital role in the star formation process by providing support against collapse, giving rise to distinct statistical properties which are used in star formation models, and by providing high density, filamentary structures where star-forming cores are preferentially located \citep{Scalo1998,Ferriere2001,MacLow2004,Kainulainen2009,Arzoumanian2011,Federrath2012,Schneider2012,Andre2014,Konstandin2016,Konyves2015,Federrath2016,Hacar2018,Mocz2018,Arzoumanian2019}. Understanding the structure, kinematics and the statistics (density and velocity dispersions, for example) of the MCs has therefore been of interest. The aim of this study is to extend upon our recent effort in \cite{Beattie2019}, herein called BFK19, to tie the fractal dimension, $\D$, to the physical properties of the MCs, e.g. to cloud length scales, and to the relations on observational data, e.g. between 2D cloud projections and 3D cloud data.

In this study we present a new method for calculating the turbulent Mach number of the clouds, based purely upon two-dimensional (2D) projected cloud position-position (PP) data, i.e., the column density, $\Sigma$, which can be obtained by molecular lines, dust emission, or dust extinction observations. We also provide the first tests of the fractal dimension methods introduced by BFK19 using dust column density maps obtained from $Herschel$ flux maps of the Polaris Flare. First, we will discuss the importance of the turbulent Mach number, and the diverse role it plays in star formation.

\subsection{The Cloud Density and Turbulent Mach Number}\label{sec:intro_Mach}
The turbulent Mach number, $\M$, is a key ingredient for numerous star formation models \citep{Krumholz2005,Federrath2010,Hennebelle2011,Federrath2012,Federrath2013,Konstandin2016}. We make the distinction between the scale-dependent turbulent Mach number,
\begin{equation} \label{eq:Mach}
\Ml(\ell) = \sigma_v(\ell)/c_s,
\end{equation}
where $\sigma_v(\ell)$ is the velocity dispersion of the cloud on length scale $\ell$, and $c_s$ is the sound speed, and the root-mean-squared (rms) Mach number,
\begin{equation} \label{eq:rmsMach}
\M \approx \sigma_v(L)/c_s = \Ml(L),
\end{equation}
where $L$ is the cloud diameter, which corresponds to the outer scale of turbulence in our study \citep{Federrath2013}. For an isothermal cloud with purely turbulent dynamics, $\M$ sets the width of the log-normal cloud density distribution,
\begin{equation}
\sigma^2_s = \ln \left(1+b^2 \M^2 \right),
\end{equation}
where the $s$ subscript denotes the variance of the normalised cloud density, $s = \ln(\rho/\rho_0)$, where $\rho_0$ is the mean density of the cloud and $b$ is the turbulent forcing parameter \citep{Padoan1997,Passot1998,Kritsuk2007,Federrath2008,Federrath2010,Konstandin2012b}. The dispersion has been studied extensively and there have been many modifications to account for 2D projections of the 3D cloud \citep{Burkhart2012}, thermal and magnetic pressures \citep{Padoan1997,Passot1998,Federrath2008,Price2011,Molina2012,Gazol2013}, and non-isothermal \citep{Nolan2015} and polytropic gases \citep{Passot1998,Li2002,Federrath2015}. Calculating the density dispersion is important for star formation models that predict the star formation rate (SFR) directly from integrating the density and free-fall time weighted cloud density distribution to determine the mass fraction of the cloud that could collapse into new stars \citep{Krumholz2005,Padoan2011,Hennebelle2011,Federrath2012,Kainulainen2014}. Beyond the density distribution, the turbulent Mach number may also play an important role in the distribution of cloud filament widths.

\subsection{Filaments and the Turbulent Mach Number}\label{sec:intro_fil}
Filament structures have been observed in star-forming and quiescent clouds, and may play an important role in star formation, since star clusters and star-forming cores have been found to be preferentially located in them \citep{Andre2010,Menschikov2010,Schneider2012,Andre2014,Padoan2014,Arzoumanian2015,Federrath2016}. Interstellar filament widths are distributed with a peak at $\sim 0.1$pc, which seems to be a universal feature of filaments and has been found in both observations and simulations of star-forming clouds \citep{Arzoumanian2011,Juvela2012,Palmeirim2013,Andre2014,Smith2014,Benedettini2015,Kirk2015,Federrath2016,Federrath2016b,Smith2016,Arzoumanian2019}. The standard deviation of the filament width distribution is thought to be associated with the sonic scale, $\ell_s$, in the clouds \citep{Federrath2018}. The sonic scale marks the transition between supersonic and subsonic velocity dispersions, and is theorised to be at length scale,
\begin{equation}
\ell_s = L \left[\frac{1}{\M} (1 + \beta)^{1/2} \right]^{1/p},
\end{equation}
in the cloud, where $p \approx 1/2$ has been measured using both Galactic cloud observations and simulations, and $\beta = p_{\text{thermal}}/p_{\text{magnetic}}$ is the ratio between the thermal and magnetic pressures at the cloud diameter scale $L$  \citep{Larson1981,Solomon1987,Ossenkopf2002,Heyer2004,Kritsuk2007,Schmidt2009,Federrath2010,RomanDuval2011,Federrath2012,Federrath2016,Federrath2018}. Being able to measure the turbulent Mach number is thus essential for testing theories about the filament width distribution. The turbulent driving that the clouds undergo may lead to the formation of filaments through interacting planar shocks \citep{Federrath2016,Tokuda2018}. However these are not the only structures that are formed and the densities of turbulent clouds have been shown to respect a fractal geometry.

\subsection{Fractal Cloud Structures}
Molecular clouds have a highly complex structure which includes sheets, filaments and dense cores. Observations of MCs through CO lines, dust emission as well as dust extinction show that they are organised into self-similar fractal structures, with substructures of clouds being continuously resolved, even at the highest spatial resolution achievable \citep{Falgarone1996,Stutzki1998,Chappell2001,Kauffmann2010,Schneider2013,Kainulainen2014,Rathborne2015}. There is a strong agreement between simulations and observations that the three-dimensional (3D) fractal dimension, i.e. the fractal dimension of the position-position-position (PPP) data of turbulent MCs falls between $\approx$ 2 and 3 (\citealt{Scalo1990,Elmegreen1996,Sanchez2005,Kowal2007,Federrath2009,Roman-Duval2010,Donovan-Meyer2013,Konstandin2016}; BFK19). However, where it falls between 2 and 3 depends upon the type of turbulent driving \citep{Federrath2009}, the rms $\M$ (\citealt{Konstandin2016}; BFK19), the length scales in the clouds, and the amount of shocks and filamentary structures in the cloud (BFK19). BFK19 also found that the fractal dimension is significantly higher in the column density map ($\approx 1.6$ in the high $\M$ limit) compared to 2D density slices ($\approx 1$ in the high $\M$ limit). In this study we show how by expanding upon the methods outlined in BFK19 one can utilise the fractal structure of the column density from the cloud, specifically the mass-length scaling, to measure $\M$, which is demonstrably an important quantity in star formation.\\

This study is organised into the following sections: In \S \ref{sec:sims} we discuss the six cloud simulations that we use to construct our new Mach number - fractal dimension $(\Ml - \D)$ relation. In \S \ref{sec:MLmethod} we summarise the fractal dimension method introduced by BFK19, including the key results. Next, in \S \ref{sec:MachNumberRelation} we derive the new relation, and discuss the applications and limitations. Then in \S \ref{sec:Application} we apply it to two quiescent subregions of the Polaris Flare to calculate the $\M$ based purely on the fractal geometry of the column density. We compare this with previous estimates of the $\M$ calculated in \cite{Schneider2013}. Finally, in \S \ref{sec:conclusion} we summarise our key findings.

\begin{figure*}
\centering
\includegraphics[width=\linewidth]{Figure1.jpg}
\caption{The column densities, $\Sigma$, for the six simulations that we use to study the dependence of the turbulent Mach number, $\M$, on the fractal dimension, $\D$. Indicated in white is the root-mean-squared (rms) $\M$ of the simulation, ranging from the transonic $\M = 1$ clouds to the highly supersonic $\M = 100$ clouds. All column densities are shown at $t = 2\,T$, in the regime of fully-developed turbulence. The densities are shown in units of log average column density, $\Sigma_0$. An animation of the time evolution for the 2D projections (and 2D slices) is available in the online version of \protect\cite{Beattie2019}.}
\label{fig:Figure1}
\end{figure*}

\section{Turbulent Molecular Cloud Models} \label{sec:sims}
In this study we use six purely hydrodynamical simulations of quiescent (non-star-forming) molecular clouds, with no self-gravitation, to construct our $\M$ -- $\D$ relation. The parameter set of the molecular cloud models is listed in Table \ref{tb:simtab}. For each of the simulations we solve the compressible Euler equations,
\begin{align}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) &= 0, \\
\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} &= - \frac{1}{\rho} \nabla P + \mathbf{F},
\end{align}
where $\rho$ is the density, $\mathbf{v}$ the velocity, $P$ the pressure, following an isothermal equation of state, $P = c_s^2 \rho$, where $c_s$ is the speed of sound, and $\mathbf{F}$ is a Ornstein-Uhlenbeck (OU) forcing function that drives the turbulence through a mixture of solenoidal $(\nabla \cdot \textbf{F} = 0)$ and compressive $(\nabla \times \mathbf{F} = 0)$ modes. We choose a natural mixture of the two modes, $b \sim 0.4$. For further details we refer the reader to \cite{Federrath2010}, \cite{Federrath2018} and BFK19.

\begin{table}
\caption{Simulation parameters.}
\centering
\begin{tabular}{r@{}lccr@{}l}
\hline
\hline
\multicolumn{2}{c}{$\M$}& Native Simulation & Number of & \multicolumn{2}{c}{Time}\\
\multicolumn{2}{c}{$(\pm1\sigma)$} & Grid Resolution & Time Slices & \multicolumn{2}{c}{Interval} \\
\hline
1.01\, & $ \pm$ 0.05 & $1024^3$ & 71 & 2 $\leq t/T$&$\leq$ 9 \\
4.1\, & $ \pm$ 0.2 & $10048^3$ & 71 & 2 $\leq t/T$&$\leq$ 9  \\
10.2\, & $ \pm$ 0.5 & $1024^3$ & 71 & 2 $\leq t/T$&$\leq$ 9\\
20\, & $ \pm$ 1 & $1024^3$ & 71 & 2 $\leq t/T$&$\leq$ 9 \\
40\, & $ \pm$ 2 & $1024^3$ & 71 & 2 $\leq t/T$&$\leq$ 9 \\
100\, & $ \pm$ 5  & $1024^3$ & 71 & 2 $\leq t/T$&$\leq$ 9 \\
\hline
\hline
\end{tabular} \\
\begin{tablenotes}
\item{\textit{\textbf{Notes:}} Column (1): the rms turbulent Mach number of the simulation $\pm$ the 1$\sigma$ temporal fluctuations. Column (2): the native 3D grid resolution of each of the simulations. Column (3): the number of time slices that we use for temporal averaging. Column (4): the time interval in units of large-scale turnover times.}
\end{tablenotes}
\label{tb:simtab}
\end{table}

\begin{figure*}
\centering
\includegraphics[width=\linewidth]{Figure2.png}
\caption{Fractal dimension as a function of length scale, for each simulation. The rms Mach numbers for each simulation are indicated in the top panels. We see low $\M$ flows associated with high fractal dimensions and clouds driven at high $\M$ to low fractal dimensions. The plot suggests that there is a smooth transition between the high and low $\M$ limits, corresponding to a transition between space-filling and shock-dominated geometries in the cloud.}
\label{fig:FDcurves}
\end{figure*}

We run the six simulations with Mach numbers $\M=1,4,10,20,40$ and $100$, over seven turnover times $(7\,T)$ in the regime of fully-developed supersonic turbulence, established after $t \geq 2\,T$ \citep{Federrath2009,Price2010}. This gives us a wide set of $\M$ values to construct the $\Ml-\D$ relation, encompassing transonic, slightly compressible flows, all the way to highly supersonic, and highly compressible flows that are saturated with shocks \citep{Federrath2013}. We solve the Euler equations in a cube with periodic boundaries (for more details on the size of the grids we refer to BFK19) and use the same initial conditions for all simulations: a homogeneous medium at rest with $\rho(x,y,z,t=0) = \text{const.}$, $\mathbf{v}(x,y,z,t=0)=0$. The same random seed for the OU forcing function is used in all simulations, hence the only difference between them is the rms Mach number.

In this study we utilise the column density data, $\Sigma$, integrated along the $z$-axis. Figure \ref{fig:Figure1} shows a single snapshot in time, $t \approx 2 \, T$ of the  column density in each simulation. In the absence of magnetic fields our turbulence simulations are isotropic in a statistical sense, i.e. when averaged over time, e.g., \cite{Federrath2009,Federrath2010}. This lets us perform our study only on the $xy$ projections (column densities), whilst still being representative 2D projections through any viewing angle.

\begin{figure*}
\centering
\includegraphics[width=\linewidth]{Figure3.png}
\caption{ \textbf{Large panel:} The turbulent Mach number as a function of fractal dimension, fit on two decades of Mach number data ($\M = 1-100$). The different colours correspond to each of the six simulations, with rms Mach number indicated in the legend. The fit (Equation \ref{eq:M(D)}) is shown in black, with 1$\sigma$ uncertainties for the fit shown as dashed lines. \textbf{Small (inset) panel:} The fractal dimension as a function of length scale for the column density. This is a key result from \protect\cite{Beattie2019}, where it is discussed in detail. The black line is the fit shown in Equation \ref{eq:D(l)}. We use the same data, but transform the length scales into Mach numbers using a power-law scaling based on supersonic turbulence theory, $\ell \sim \Ml^2$, to create the $\Ml$ -- $\D$ relation \protect\citep{Burgers1948,Konstandin2012,Federrath2013}. }
\label{fig:MachFDRelation}
\end{figure*}

\section{Fractal Dimension Curves} \label{sec:MLmethod}
We follow the mass-length fractal dimension method outlined and discussed with detail in BFK19. This method allows us to calculate a mass-length $\D$ on each length scale in the cloud. It is important that we are able to access $\D$ on each $\ell$, since our aim is to relate $\D$ with $\Ml$ through $\ell$. We provide a summary of the method, and the key results below. Please note in BFK19 we use $\D_{\text{p}}$ to indicate the 2D projected (column density) fractal dimension but in this study we use $\D$ for simplicity.

\begin{table}
\caption{Fitting parameter values for $\D(\ell/L)$, shown in Equation \ref{eq:D(l)} and $\D(\Ml)$, shown in Equation \ref{eq:D(M)}.}
\centering
\begin{tabular}{cccccc}
\hline
\hline
$\Dmin \pm 1\sigma$ & $\Dmax$ & $\beta_0 \pm 1\sigma$ & $\beta_1 \pm 1\sigma$ & $\M \pm 1\sigma$ \\
\hline
$1.55 \pm 0.13$ & 2 & $0.46 \pm 0.08$ & $0.56 \pm 0.08$  & $4.1 \pm 0.2$ \\
\hline
\hline
\end{tabular} \\
\begin{tablenotes}
\item{\textit{\textbf{Notes:}} $\beta_0$, $\beta_1$ and $\Dmin$ are calculated in \cite{Beattie2019} using weighted non-linear regression. Column (1): $\Dmin$ is the minimum fractal dimension of the column density. Column (2): $\Dmax$ is assumed to be 2 for the maximum fractal dimension of the column density. This corresponds to completely space-filling flows on the 2D plane. Column (3): $\beta_0$ is a fitting parameter that corresponds to the translation of the complimentary error function over the $\ell/L$ axis. Column (4): $\beta_1$ is a fitting parameter that corresponds to the rate in which the complimentary error function changes between the high $\D$ and low $\D$ states. Column (5): The rms Mach number, $\M$ is measured by averaging over all rms Mach numbers from $2 \leq t/T \leq 9$ in the $\M=4$ simulation. The $1\sigma$ uncertainty is the standard deviation of the averaging process. We use $\M=4$ because the curves are in the Mach 4 simulation frame of reference, which is discussed in detail in \cite{Beattie2019}.}
\end{tablenotes}
\label{tb:parms1}
\end{table}

\subsection{Method Summary} \label{sec:methodsummary}
There exists a power-law scaling between the mass and length scales (size) in real MCs \citep{Larson1981,Myers1983,Falgarone1996,Roman-Duval2010,Donovan-Meyer2013}. To calculate the mass-length dimension we utilise this power-law scaling, where the mass, $m$ is given by
\begin{equation}\label{eq:ML}
m(\ell/L) \sim \left(\ell/L\right)^{\D},
\end{equation}
where we use the dimensionless $\ell/L$ for the length scales in the cloud, and $\D$ for the scaling exponent, i.e. the fractal dimension for the mass-length scaling relation. The constant of proportionality is $M$, the total mass, since when $\ell = L$, $m(1) \sim 1^{\D} = M$, in the unit system used here. We explore how the fractal dimension changes with spatial scale, i.e. we examine the $\ell/L$ dependence of the scaling exponent, $\D$,
\begin{equation}\label{eq:MLS}
m(\ell/L) = M (\ell/L)^{\mathcal{D}(\ell/L)}.
\end{equation}
We use Equation \ref{eq:MLS} to then define the fractal dimension at length scale $\ell_i/L$,
\begin{equation} \label{eq:MLS2}
    \mathcal{D}(\ell_i/L) = \frac{\log(m_i/M)}{\log(\ell_i/L)},
\end{equation}
where the $\ell_i/L$ length scale can be interpreted as all nested length scales up to the length scale $\ell_i/L$, i.e. $\ell_i / L = \left\{ \ell_0/L, \ell_1 /L , \hdots, \ell_i / L \right\}$, where $\ell_0/L < \ell_1/L < \hdots < \ell_i/L$ and $\ell_0/L$ is the smallest possible length scale in the cloud, and $m_i$ the corresponding mass at the $\ell_i / L$ scale, i.e., the total mass of the cloud on all scales less than and including $\ell_i / L$. This treats the cloud like a nested, hierarchical set of density objects, each with its own $\mathcal{D}$ and lets us probe how self-similar the cloud is over all length scales. If it is self-similar over a set of $\ell/L$, for example, the fractal dimension will not change as a function of $\ell/L$ over this region. This method also can be easily extended to explore the scale-dependent density structure of the cloud,
\begin{align}
    \frac{m}{(\ell/L)^2} &\sim (\ell/L)^{\D -2}, \\
    \Sigma &\sim (\ell/L)^{\D - 2},
\end{align}
to access the scaling in the column density, $\Sigma$.

To calculate the $\D$ described above we need to calculate the mass as a function of length scale. We do this by performing the following steps on the $\Sigma$ data from each simulation, following exactly the method outlined BFK19:

\begin{enumerate}
    \item Identify the coordinates of the maximum column density pixel $\Sigma_{\text{max}}$, \\[0.5em]
    \item expand a $(\ell/L)\times(\ell/L)$ square region centred on $\Sigma_{\text{max}}$, creating our length scale hierarchy, \\[0.5em]
    \item calculate the mass $m$, within each of the $(\ell/L)\times(\ell/L)$ squares, \\[0.5em]
    \item use the relation shown in Equation \ref{eq:MLS2} to determine the fractal dimension, $\D$, as a function of $\ell_i/L$, always using all length scales below $\ell_i/L$ to calculate $\D$ on $\ell_i/L$. \\[0.5em]
\end{enumerate}

We apply the four steps above on each of the 71 time slices in the interval $2 \leq t/T \leq 9$, the statistically fully-developed turbulence regime, averaging over them to construct a single curve for each simulation with 1$\sigma$ uncertainties.

\subsection{Key Results from the $\D(\ell)$ Curves}
In BFK19 we show that the fractal dimension curves from each simulation (shown in Figure \ref{fig:FDcurves}) can be combined together into the same common reference frame to create a composite fractal dimension curve. This lets us map the fractal dimension of $\Sigma$ over seven orders of magnitude in spatial scales, encompassing clouds undergoing subsonic to highly supersonic turbulent dynamics. After combining the curves we find that a complimentary error function is a good fit for $\D(\ell/L)$, which models a smooth transition between space-filling clouds ($\D=2$ in the 2D projection) and clouds saturated with planar shocks ($\D=1.55\pm 0.13$, which is a key result from BFK19). A simple power-law relation is not sufficient, because at both low and high rms Mach flows the fractal dimension curves begin to flatten out as they approach the high and low limits. The empirical fit is,
\begin{equation}\label{eq:D(l)}
\D(\ell/L) = \frac{\Dmax - \Dmin}{2}\erfc \Big( \beta_1\log \left( \ell/L \right) + \beta_0 \Big) + \Dmin,
\end{equation}
where $\erfc$ is the complimentary error function, $\Dmin$ and $\Dmax$ are the minimum and maximum fractal dimension, respectively, and $\beta_1$ and $\beta_0$ are fitting parameters, determined using nonlinear least squares and tabulated in Table \ref{tb:parms1}. This fit encodes the limits $\lim_{(\ell/L) \rightarrow  1 } \D(\ell/L) = \Dmin$ and $\lim_{(\ell/L) \rightarrow  0} \D(\ell/L) = \Dmax$ that we find in the data, and encompasses the smooth transition that we find between them. The composite curve data, along with the complimentary error function fit for $\D(\ell/L)$ are shown in the sub-panel of Figure \ref{fig:MachFDRelation}. The plot shows that the $\D$ of the column density is bounded between $\Dmax = 2$ and $\Dmin = 1.55 \pm 0.13$, where the former is assumed, and the latter is measured through the fitting process. Next, we use this curve to construct the $\Ml - \D$ relation.

\begin{figure*}
\centering
\includegraphics[width=\linewidth]{Figure5.png}
\caption{The dust column density maps of the `saxophone' and `quiet' subregions of the Polaris Flare shown in $\log_{10}$ [cm$^{-22}$]. The maximum column density pixels, $\Sigma_\text{max}$ are indicated on the maps. These pixels were used to initialise the $\ell\times\ell$ expansion for the construction of $\D(\ell)$. The fractal dimensions of the regions are calculated to be $\D = 1.60 \pm 0.04$ and $\D = 1.76 \pm 0.05$, for `saxophone' and `quiet', respectively. This is consistent with \protect\cite{Beattie2019} where it is argued that low $\D$ values are associated with high $\M$ clouds, and high $\D$ values are associated with low $\M$ clouds.}
\label{fig:fig5}
\end{figure*}

\section{The $\Ml$ -- $\D$ Relation} \label{sec:MachNumberRelation}
After establishing the length scale dependence of the fractal dimension we may immediately ask how then does the fractal dimension change with the velocity dispersion of the cloud, since the velocity dispersion also depends upon length scale, as indicated in Equation \ref{eq:Mach}. This link is made available to us by approximating $\Ml(\ell)$ using scaling relations from models of supersonic turbulence.

\subsection{Constructing the Relation}
Using the second-order structure function, $\SF_2(\ell/L) = \left\langle | \mathbf{v}(\mathbf{r}) - \mathbf{v}(\mathbf{r} + \ell/L)|^2 \right\rangle_{\mathbf{r}},$ where $\mathbf{v}(\mathbf{r})$ is the velocity of the cloud at position $\mathbf{r}$, and the operator $\left\langle \hdots \right\rangle_{\mathbf{r}}$ is an average over a large ensemble of spatial positions, one can construct the turbulent Mach number at length scale $\ell/L$ (using the definition of the second-order structure function and Equation 4.3 from \citealt{Konstandin2012}). This construction of the Mach number follows a power-law of the form,
\begin{equation}
\Ml(\ell/L) = \sqrt{\SF_2(\ell/L)/(2c_s^2)} \sim (\ell/L)^p,
\end{equation}
where $p \approx 1/2$ for supersonic turbulence and $p \approx 1/3$ for subsonic turbulence \citep{Kolmogorov1941,Burgers1948,Kritsuk2007,Schmidt2009,Konstandin2012,Federrath2013,Federrath2018}.
We use the $p = 1/2$ case, $\ell/L \sim \Ml(\ell/L)^{2}$, to transform all length scales into Mach numbers. We set $\ell = L$ to find the constant of proportionality, $\M^{-2}$, i.e. on large scales in the cloud $\Ml$ is $\approx \M$ \citep{Federrath2013}. Hence the transformation is

\begin{equation} \label{eq:Mlmap}
\ell/L = [\Ml(\ell/L)/\M]^2.
\end{equation}

We apply this transformation to our $\ell/L$ composite data shown in the sub-panel of Figure \ref{fig:MachFDRelation}. This provides us with an estimate for $\mathcal{D}(\Ml)$,

\begin{equation}\label{eq:D(M)}
\D(\Ml) = \frac{\Dmax - \Dmin}{2}\erfc \Bigg[ 2 \beta_1 \log \left( \frac{\Ml}{\M} \right) + \beta_0 \Bigg] + \Dmin.
\end{equation}

We invert the equation to obtain $\Ml(\mathcal{D})$, which means that from measurements of the fractal dimension one can infer the scale-dependent turbulent Mach number,

\begin{equation}\label{eq:M(D)}
\log\Ml(\D) = \xi_1\left(\erfc^{-1}\left[ \frac{\D-\Dmin}{\Omega} \right] + \xi_2 \right),
\end{equation}
where
$$\Omega = \frac{\Dmax-\Dmin}{2},$$
and
$$\xi_1 = (2\beta_1)^{-1} \hspace{2mm} \text{and} \hspace{2mm} \xi_2 = \xi_1^{-1}\log\M - \beta_0.$$
This corresponds to the turbulent Mach number on the length scale that $\D$ was measured, since $\Ml(\D(\ell))$. The values of the estimated and derived parameters for this fit are shown in Table \ref{tb:parms2}.

\begin{table}
\centering
\caption{Derived parameter values for $\Ml(\D)$, shown in Equation \ref{eq:M(D)}.}
\begin{tabular}{cccc}
\hline
\hline
$\Dmin \pm 1\sigma$ & $\Omega \pm 1\sigma$ & $\xi_1 \pm 1\sigma$ & $\xi_2 \pm 1\sigma$ \\
\hline
$1.55 \pm 0.13$ & $0.22 \pm 0.07$ & $0.9 \pm 0.1$ & $0.2 \pm 0.2$ \\
\hline
\hline
\end{tabular} \\
\begin{tablenotes}
\item{\textit{\textbf{Notes:}} The values for $\Dmax$, $\beta_0$, $\beta_1$ and $\M$ can be found in Table \ref{tb:parms1}, which are used to derive the parameters in this table. Column (1): the same as column one from Table \ref{tb:parms1}, but included for completeness. Column (2): $\Omega = (\Dmax - \Dmin)/2$. Column (3): $\xi_1 = (2\beta_1)^{-1}$. Column (4): $\xi_2 =  \xi^{-1}_1 \log \M - \beta_0$.}
\end{tablenotes}
\label{tb:parms2}
\end{table}

\subsection{$\Ml$ -- $\D$ Relation Results}
In the main plot of Figure \ref{fig:MachFDRelation} we show $\Ml$ as a function of $\D$, derived from the column density data. The black line shows Equation \ref{eq:M(D)}, which is the equation we will use to convert fractal dimensions into Mach numbers. There are three main limitations to this method. The first is that the inverse complimentary error function has steep tails. This means that for low and high $\M$ the relation is extremely sensitive to small changes in $\D$. The second is that for high $\M$ the temporal fluctuations of $\D$ become significant, spanning over $> 0.2 \D$. This means that our relation will perform best at measuring regions of MCs with $10 \gtrsim \M \gtrsim 0.1$. Finally, since in our construction of the relation we use clouds driven by compressible, supersonic, isothermal and isotropic turbulence, with a natural mixture between solenoidal and compressive modes $(b \sim 0.4)$, the fit shown in Figure \ref{fig:MachFDRelation} may only work well, without modification, on quiescent clouds, without significant deviation from natural mixing, since $\D$ is sensitive to changes in driving \cite{Federrath2009}. Acknowledging these limitations, we now test the performance of the relation on $Herschel$ observations of subregions from the Polaris Flare cloud.

\section{Application on quiescent subregions in the Polaris Flare Cloud} \label{sec:Application}
The Polaris Flare is a high Galactic latitude cloud which is located at a distance $\leq 150$pc \citep{Falgarone1998,Miville2010,Schlafly2014}. It has weak, but significant CO emissions in regions of higher hydrogen column density \citep{Falgarone1998,Meyerdierks1996,Miville2010}. There is no active star-formation in Polaris, only 5 starless cores were detected \citep{Andre2010,Thompson2010} that are most likely not gravitationally bound. Polaris is thus a perfect candidate for testing our new relation, which was calibrated upon quiescent cloud simulations.

\subsection{Observational Data}
The Polaris region was observed as part of the $Herschel$ Gould Belt survey (HGBS, \citealt{Andre2010}) using the PACS and SPIRE instruments on-board $Herschel$. For all observational details, we refer to \citep{Andre2010,Menschikov2010,Thompson2010,Miville2010}. We employ publicly available level 3 data products produced with HIPE13 (Herschel Interactive Processing Environment) from the $Herschel$ archive. The angular resolution of the maps is 11.7$''$, 18.2$''$, 24.9$''$, and 36.3$''$ for 160 $\mu$m (PACS) and 250, 350, and 500 $\mu$m (SPIRE), respectively. For an absolute calibration of the maps (included in the SPIRE level 3 data), the Planck High Frequency Instrument (HFI) observations were used for the HIPE-internal zero-point correction task that calculates the absolute offsets, based on cross-calibration with HFI-545 and HFI-857 maps, including colour-correcting HFI to SPIRE wavebands, assuming a grey-body function with fixed beta. For the PACS 160 $\mu$m map, we obtained the zero-point correction following the procedure outlined in \cite{Bernard2010}.
Column density and temperature maps were then produced at an angular resolution of 18$''$, following the procedure outlined in \cite{Palmeirim2013} that employs a multi-scale decomposition of the imaging data and assumes a constant line-of-sight temperature. We performed a pixel-by-pixel SED (Spectral Energy Distribution) fit from 160 to 250$\mu$m, using a dust opacity law $\kappa_0 = 0.1 \times (\nu/1000GHz)^{\,\beta}$ cm$^2$ g$^{-1}$ with $\beta$ = 2 and assuming a gas-to-dust ratio of 100. This dust opacity law is commonly adopted in other HGBS publications and we refer to \cite{Andre2010} and \cite{Konyves2015} for further details. We estimate that the final uncertainties of the column density map are around 20--30 \%. \\
The resulting large, high-angular resolution (18$''$) hydrogen column density map traces structures between 0.01 to 8 pc \citep{Miville2010,Schneider2013}. For our study, we cut out the two subregions that have previously been used to investigate the link between the probability distribution function of column density (N-PDF) and $\M$ in Polaris \citep{Schneider2013}, the `saxophone' and `quiet' subregions, shown in Figure \ref{fig:fig5}. \cite{Schneider2013} made estimates of the turbulent Mach number for the two regions based on the CO\footnote{The CO data stem from $^{12}$CO 2$\to$1 and $^{13}$CO 1$\to$0 observations from the KOSMA and FCRAO telescopes  published in \cite{Bensch2003}.} velocity dispersion. The authors
calculated the $\M$ for each subregion using,
\begin{equation} \label{eq:S.Mach}
    \M = \left(\sqrt{3}\text{FWHM} \right) \Big/ \left(c_s \sqrt{8\ln 2} \right)
\end{equation}
where FWHM [km s$^{-1}$] is the full width at half maximum of the CO molecular line data, and under the LTE assumption the sound speed is $c_s \approx 0.188\sqrt{\text{T}_{\text{ex}}/10\text{K} }$, where $\text{T}_{\text{ex}}$ is the excitation temperature, $\text{T}_{\text{ex}} = 5.53\left[ \ln(5.33/\text{T}_{\text{CMB}}) + 1 \right]^{-1}$ and $\text{T}_{\text{CMB}} = 2.728$K is the temperature of the cosmic microwave background. The values for these estimates of $\M$ are shown in column (2) of Table \ref{tb:MachMeasurements}.

\subsection{Application and considerations}
Using the method summarised in \S \ref{sec:methodsummary} we construct the fractal dimension curves for each of the two regions and then use Equation \ref{eq:M(D)} to calculate the Mach number\footnote{Our fractal dimension implementation on the Polaris Flare cloud is available here: \url{https://github.com/AstroJames/FractalGeometryofPolaris}}. The application differs in two ways compared to BFK19. (1) Here we use terminating boundary conditions, while on the simulation data we use periodic boundaries. Since the observational data is not periodic we terminate the $\ell\times\ell$ expansion on the boundaries. (2) We calculate the monofractal mass-length dimension instead of our length-dependent curves. Our method reduces exactly to the monofractal $\D$ when $\D(\ell = L)$, where $L$ is the largest scale of the $\ell\times\ell$ expansion. This means that we use all nested scales $\ell \leq L$ to calculate $\D$, as described in \S \ref{sec:methodsummary}. We do this because it allows us to calculate $\M$, the turbulent Mach number on the cloud scale, which is useful for constraining the 3D density PDF, calculating the star formation rate and estimating the sonic scale in the cloud (as discussed in \S \ref{sec:intro_Mach} and \ref{sec:intro_fil}, respectively). It also allows for comparison with previous estimates made for $\M$ in the `quiet' and `saxophone' subregions of the Polaris Flare molecular cloud using Equation \ref{eq:S.Mach}.

\subsection{$\M(\D)$ Results}
In Figure \ref{fig:fig4} we show our $\M$ estimates for the two subregions of Polaris, indicated in red, and compare them with the $\M$ estimated by \cite{Schneider2013}, shown as blue-dashed regions. The $\M$ calculated by \cite{Schneider2013} and the $\M$ calculated in this study are shown in Table \ref{tb:MachMeasurements}. We find $\M = 2\pm 1$ and $\M = 10^{+55}_{-4}$, whereas \cite{Schneider2013} finds $\M 3 \pm 1$ and $\M = 7 \pm 3$, for the `quiet' and `saxophone' subregions of Polaris, respectively. These estimates are consistent to within 1$\sigma$. However the Mach number measurements are very sensitive to small changes in $\D$, especially for $\M \gtrsim 10$, where the relation becomes extremely steep. This translates into large, and not necessarily symmetric uncertainties, as shown for the `saxophone' subregion, in Table \ref{tb:MachMeasurements}. To understand why there may be differences between the $\M(\D)$ and previous estimates based on the CO velocity dispersion we turn to the calculation of $\D$.

To estimate $\M$ we first calculate $\D$, which is shown in column (3) of Table \ref{tb:MachMeasurements}. We find `quiet', the lower $\M$ region, has a $\D = 1.76 \pm 0.05$, and `saxophone', the higher $\M$ region, $\D = 1.60 \pm 0.04$. This is consistent with BFK19, who argues that for higher $\M$ flows we should expect lower $\D$, corresponding to the introduction of compressive shocks into a diffuse cloud with increasing $\M$, and previous studies have calculated $\D = 1.4 - 1.8$ for column densities \citep{Elmegreen1996,Elmegreen2004,Sanchez2005,Rathborne2015}. This suggests that the $\D$ values we calculate are reasonable, however, the `saxophone' region has a clear filamentary structure (see the high-density filament feature on the left column density map in Figure \ref{fig:fig5}), which has a density-length scaling relation $\rho \sim \ell^{-2}$, or $\Sigma \sim \ell^{-1}$, and will act to reduce $\D$ in the vicinity of $\Sigma_{\text{max}}$ \citep{Schneider2013,Federrath2016,Andre2017}. This may account for the slightly higher $\M$ that we estimate for `saxophone', compared to the Mach estimate in \cite{Schneider2013}. For the `quiet' subregion we slightly underestimate $\M$. This is because the column density (see the right column density map in Figure \ref{fig:fig5}) is diffuse, and lacks the shock structures that we see introduced between the $\M = 1$ and $\M = 4$ simulations in Figure \ref{fig:Figure1}. The different column density geometry in the `quiet' cloud may be due to a deviation away from the natural mixing of driving modes, towards stronger, compressive driving which we do not currently include in our relation, and which can change the fractal dimension up to $\sim 0.1\D$ for the mass-length method \citep{Federrath2009}.

\begin{table}
\caption{Mach number values for both Polaris Flare subregions.}
\centering
\begin{tabular}{cccc}
\hline
\hline
Subregion  & $\M$ (based on Eq. \ref{eq:S.Mach}) & $\D \pm 1\sigma$  & $\M(\D) \pm 1\sigma$  \\
\hline
Saxophone       & $ 7 \pm 3$  & $1.60 \pm 0.04$   & $10^{+55}_{-4}$ \\
Quiet           & $3 \pm 1$  & $1.76 \pm 0.05$   & $2 \pm 1$ \\
\hline
\hline
\end{tabular} \\
\begin{tablenotes}
\item{\textit{\textbf{Notes:}} Column (1): The subregion in the Polaris Flare cloud. Column (2): the estimated $\M$ from \cite{Schneider2013}, which is stated to have an error $\sim$30-40\%. Column (3): $\D$ calculated from our mass-length method. We take the $\D$ on the largest scale calculated, which reduces to the regular monofractal mass-length dimension. Column (4): The Mach number calculated from Equation \ref{eq:M(D)}, with 1$\sigma$ uncertainties propagated from column (3).}
\end{tablenotes}
\label{tb:MachMeasurements}
\end{table}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{Figure4.png}
\caption{The $\Ml$ -- $\D$ relation shown in black, with $\M$ measurements of the quiet and saxophone regions in the Polaris Flare cloud shown as red markers. The light-blue bands are the propagated $1\sigma$ uncertainties from the $\D$ measurements. The $\M_{\text{sax}} \sim 7$ and $\M_{\text{qui}} \sim 3$ regions, encompassed by the blue, dashed lines indicate the measured $\M$ using the CO velocity dispersion of the regions made by \protect\cite{Schneider2013}. We calculate $\M \sim 10$ and $\M \sim 2$ for the two regions, respectively, which is consistent with \protect\cite{Schneider2013} to within 1$\sigma$. The values of all $\M$ and $\D$ measurements can be found in Table \ref{tb:MachMeasurements}. }
\label{fig:fig4}
\end{figure}

\section{Summary and Key Findings}\label{sec:conclusion}
In this study we construct a new empirical relation for the scale-dependent three-dimensional (3D) Mach number, $\Ml$, and the fractal dimension, $\D$, of the column density for turbulent clouds. We use the mass-length fractal dimension method introduced in \cite{Beattie2019} (BFK19) on six hydrodynamical cloud simulations, with root-mean-squared (rms) Mach number, $\M$, varying from $\M = 1$ to $100$. We apply the method on the column densities, with an example of the densities shown Figure \ref{fig:Figure1}, to construct $\D$ as a function of length scale, $\ell$. We then transform the cloud length scales to $\Ml$ using the scaling relation $\ell \sim \Ml^2$, for supersonic turbulence \citep{Burgers1948,Federrath2013}. Using this data we are able to construct $\D(\Ml)$, and finally $\Ml(\D)$. Using $\Ml(\D)$ we construct $\M$, where $\M \sim \Ml(L)$, and $L$ is the cloud scale, for the dust column density maps of two quiescent subregions from the Polaris Flare, that are termed `quiet' and `saxophone', studied earlier in \cite{Schneider2013}. We summarise our key findings below:

\begin{itemize}
\item We propose a new empirical relation for the scale-dependent Mach number and the fractal dimension of the column density, $$\log\Ml(\D) = \xi_1\left(\erfc^{-1}\left[ \frac{\D-\Dmin}{\Omega} \right] + \xi_2 \right),$$ as defined in Equation \ref{eq:M(D)} and plotted in Figure \ref{fig:MachFDRelation}, where $\xi_1 = 0.9 \pm 0.1$, $\xi_2 = 0.2 \pm 0.2$, $\Omega = 0.22 \pm 0.07$ and $\Dmin$, the fractal dimension of the column density in the high $\M$ limit is $1.55 \pm 0.13$. This relation allows for the calculation of $\M$ for clouds in the range $10 \gtrsim \M \gtrsim 0.1$. Very large and very low $\M$ are inappropriate for the model due to the steep tails in the inverse complimentary error function. \\[0.2em]
\item We use the mass-length fractal dimension method in BFK19 to calculate the $\D$ of the `saxophone' and `quiet' subregions of the Polaris Flare, shown in Figure \ref{fig:fig5}. We find $\D = 1.60 \pm 0.04$ and $\D = 1.76 \pm 0.05$ for `saxophone' and `quiet', respectively, consistent with the thesis that higher Mach number flows reduce $\D$, by turning diffuse, space-filling structures into compressive shocks and filaments. \\[0.2em]
\item Using $\D$ we estimate $\M$ for each of the subregions. We find `quiet' has a $\M \sim 2$ and `saxophone' has a $\M \sim 10$, shown in Figure \ref{fig:fig4}. This is comparable to the estimates made in \cite{Schneider2013}, $\M \sim 3$ and $\M \sim 7$ for the two respective subregions, but based on the CO velocity dispersion. The agreement between the Mach number estimate based on the CO velocity dispersion and based on our new fractal dimension relation, Equation \ref{eq:M(D)}, is acceptable, especially considering we do not account for how different types of turbulent driving influence the cloud geometries or how the presence of large filamentary structures, that locally scale the cloud by $\Sigma \sim \ell^{-1}$, act to reduce $\D$. \\[0.2em]
\item Our results suggest that the new empirical relation between the fractal dimension of column densities and 3D turbulent Mach number is a useful tool for extracting the Mach number purely from the structure and geometry of column density data from the cloud.
\end{itemize}

\section*{Acknowledgements}\label{section:acknowledgments}
We thank the anonymous reviewer for the detailed and critical reading, which improved the study. J.~R.~B~acknowledges Tilly for the ever-present support. C.~F.~acknowledges funding provided by the Australian Research Council (Discovery Project DP170100603, and Future Fellowship FT180100495), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). R.~S.~K.~acknowledges support from the Deutsche Forschungsgemeinschaft via SFB 881, ``The Milky Way System" (sub-projects B1, B2 and B8). We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants~pr32lo, pr48pi and GCS Large-scale project~10391), the Partnership for Advanced Computing in Europe (PRACE grant pr89mu), the Australian National Computational Infrastructure (grant~ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, in the framework of the National Computational Merit Allocation Scheme and the ANU Allocation Scheme. The simulation software FLASH was in part developed by the DOE-supported Flash Centre for Computational Science at the University of Chicago. N.~S.~acknowledges support by the French ANR and the German DFG through the project "GENESIS" (ANR-16-CE92-0035-01/DFG1591/2-1).

\bibliographystyle{mnras.bst}