The relation between the turbulent Mach number and observed fractal dimensions of turbulent clouds

TL;DR

This paper introduces a new method to estimate the turbulent Mach number of molecular clouds using their 2D column density fractal dimension, enabling kinematic analysis from purely geometric data.

Contribution

The authors develop an empirical relation between fractal dimension and Mach number, validated with simulations and Herschel observations, facilitating Mach number estimation from 2D data.

Findings

The empirical relation accurately estimates Mach numbers from column density maps.

Application to Herschel data yields Mach numbers consistent with traditional velocity dispersion methods.

The method enables kinematic insights solely from geometric column density measurements.

Abstract

Supersonic turbulence is a key player in controlling the structure and star formation potential of molecular clouds (MCs). The three-dimensional (3D) turbulent Mach number, $\mathcal{M}$, allows us to predict the rate of star formation. However, determining Mach numbers in observations is challenging because it requires accurate measurements of the velocity dispersion. Moreover, observations are limited to two-dimensional (2D) projections of the MCs and velocity information can usually only be obtained for the line-of-sight component. Here we present a new method that allows us to estimate $\mathcal{M}$ from the 2D column density, $\Sigma$, by analysing the fractal dimension, $\mathcal{D}$. We do this by computing $\mathcal{D}$ for six simulations, ranging between $1$ and $100$ in $\mathcal{M}$. From this data we are able to construct an empirical relation, $\log\mathcal{M}(\mathcal{D}) = \xi_1(\text{erfc}^{-1} [(\mathcal{D}-\mathcal{D}_{\text{min}})/\Omega] + \xi_2),$ where $\text{erfc}^{-1}$ is the inverse complimentary error function, $\mathcal{D}_{\text{min}} = 1.55 \pm 0.13$ is the minimum fractal dimension of $\Sigma$, $\Omega = 0.22 \pm 0.07$, $\xi_1 = 0.9 \pm 0.1$ and $\xi_2 = 0.2 \pm 0.2$. We test the accuracy of this new relation on column density maps from $Herschel$ observations of two quiescent subregions in the Polaris Flare MC, `saxophone' and `quiet'. We measure $\mathcal{M} \sim 10$ and $\mathcal{M} \sim 2$ for the subregions, respectively, which is similar to previous estimates based on measuring the velocity dispersion from molecular line data. These results show that this new empirical relation can provide useful estimates of the cloud kinematics, solely based upon the geometry from the column density of the cloud.

Tables

Table 1. Table 1: Simulation parameters.

$ℳ$		Native Simulation	Number of	Time
$\left(\pm 1\sigma \right)$		Grid Resolution	Time Slices	Interval
1.01	$\pm$ 0.05	${1024}^3$	71	2 $\le t/T$	$\le$ 9
4.1	$\pm$ 0.2	${10048}^3$	71	2 $\le t/T$	$\le$ 9
10.2	$\pm$ 0.5	${1024}^3$	71	2 $\le t/T$	$\le$ 9
20	$\pm$ 1	${1024}^3$	71	2 $\le t/T$	$\le$ 9
40	$\pm$ 2	${1024}^3$	71	2 $\le t/T$	$\le$ 9
100	$\pm$ 5	${1024}^3$	71	2 $\le t/T$	$\le$ 9

Table 2. Table 2: Fitting parameter values for 𝒟 ⁡ ( ℓ / L ) 𝒟 ℓ 𝐿 \operatorname{\mathcal{D}}(\ell/L) , shown in Equation 12 and 𝒟 ⁡ ( ℳ ) 𝒟 ℳ \operatorname{\mathcal{D}}(\operatorname{\mathscr{M}}) , shown in Equation 15 .

${𝒟}_{\text{min}}\pm 1\sigma$	${𝒟}_{\text{max}}$	${\beta }_0\pm 1\sigma$	${\beta }_1\pm 1\sigma$	$ℳ\pm 1\sigma$
$1.55\pm 0.13$	2	$0.46\pm 0.08$	$0.56\pm 0.08$	$4.1\pm 0.2$

Table 3. Table 3: Derived parameter values for ℳ ⁡ ( 𝒟 ) ℳ 𝒟 \operatorname{\mathscr{M}}(\operatorname{\mathcal{D}}) , shown in Equation 16 .

${𝒟}_{\text{min}}\pm 1\sigma$	$\Omega \pm 1\sigma$	${\xi }_1\pm 1\sigma$	${\xi }_2\pm 1\sigma$
$1.55\pm 0.13$	$0.22\pm 0.07$	$0.9\pm 0.1$	$0.2\pm 0.2$

Table 4. Table 4: Mach number values for both Polaris Flare subregions.

Subregion	$ℳ$ (based on Eq. 17)	$𝒟\pm 1\sigma$	$ℳ\left(𝒟\right)\pm 1\sigma$
Saxophone	$7\pm 3$	$1.60\pm 0.04$	${10}_{-4}^{+55}$
Quiet	$3\pm 1$	$1.76\pm 0.05$	$2\pm 1$

Equations

$$M(\ell )={\sigma }_v(\ell )/c_s,$$

$$M\approx {\sigma }_v(L)/c_s=M(L),$$

$${\sigma }_s^2=\ln \left(1+b^2{M}^2\right),$$

$${\ell }_s=L{\left[\frac{1}{M}(1+\beta {)}^{1/2}\right]}^{1/p},$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})$$

$$\frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}\cdot \nabla )\mathbf{v}$$

$$m(\ell /L)\sim {\left(\ell /L\right)}^D,$$

$$m(\ell /L)=M(\ell /L{)}^{\mathcal{D}(\ell /L)}.$$

$$\mathcal{D}({\ell }_i/L)=\frac{\log (m_i/M)}{\log ({\ell }_i/L)},$$

$$\frac{m}{(\ell /L{)}^2}$$

$$\Sigma$$

$$\operatorname{\mathcal{D}}(\ell/L)=\frac{\operatorname{\mathcal{D}_{\text{max}}}-\operatorname{\mathcal{D}_{\text{min}}}}{2}\operatorname{erfc}\Big{(}\beta_{1}\log\left(\ell/L\right)+\beta_{0}\Big{)}+\operatorname{\mathcal{D}_{\text{min}}},$$

$$M(\ell /L)=\sqrt{{\mathrm{SF}}_2(\ell /L)/(2c_s^2)}\sim (\ell /L{)}^p,$$

$$\ell /L=[M(\ell /L)/M{]}^2.$$

$$\operatorname{\mathcal{D}}(\operatorname{\mathscr{M}})=\frac{\operatorname{\mathcal{D}_{\text{max}}}-\operatorname{\mathcal{D}_{\text{min}}}}{2}\operatorname{erfc}\Bigg{[}2\beta_{1}\log\left(\frac{\operatorname{\mathscr{M}}}{\operatorname{\mathcal{M}}}\right)+\beta_{0}\Bigg{]}+\operatorname{\mathcal{D}_{\text{min}}}.$$

$$\log M(D)={\xi }_1\left({\mathrm{erfc}}^{-1}\left[\frac{D-{\mathcal{D}}_{\text{min}}}{\Omega }\right]+{\xi }_2\right),$$

$$\Omega =\frac{{\mathcal{D}}_{\text{max}}-{\mathcal{D}}_{\text{min}}}{2},$$

$${\xi }_1=(2{\beta }_1{)}^{-1}\text{and}{\xi }_2={\xi }_1^{-1}\log M-{\beta }_0.$$

$$\operatorname{\mathcal{M}}=\left(\sqrt{3}\text{FWHM}\right)\Big{/}\left(c_{s}\sqrt{8\ln 2}\right)$$

$$\log M(D)={\xi }_1\left({\mathrm{erfc}}^{-1}\left[\frac{D-{\mathcal{D}}_{\text{min}}}{\Omega }\right]+{\xi }_2\right),$$

Full text

††pubyear: 2019††pagerange: The relation between the turbulent Mach number and observed fractal dimensions of turbulent clouds–References

1 Introduction

The dynamical evolution of molecular clouds (MCs) in the interstellar medium (ISM) is determined by supersonic, compressible turbulent flows (Larson, 1981; Solomon et al., 1987; Klessen et al., 2000; Heitsch et al., 2001; Ossenkopf & Mac Low, 2002; Elmegreen & Scalo, 2004; Heyer & Brunt, 2004; Mac Low & Klessen, 2004; Scalo & Elmegreen, 2004; Krumholz & McKee, 2005; Ballesteros-Paredes et al., 2007; McKee & Ostriker, 2007; Roman-Duval et al., 2011; Padoan et al., 2014; Federrath & Banerjee, 2015). 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 (Scalo, 1998; Ferrière, 2001; Mac Low & Klessen, 2004; Kainulainen et al., 2009; Arzoumanian et al., 2011; Federrath & Klessen, 2012; Schneider et al., 2012; André et al., 2014; Konstandin et al., 2016; Könyves et al., 2015; Federrath, 2016; Hacar et al., 2018; Mocz & Burkhart, 2018; Arzoumanian et al., 2019). 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 Beattie et al. (2019), herein called BFK19, to tie the fractal dimension, $\operatorname{\mathcal{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.

1.1 The Cloud Density and Turbulent Mach Number

The turbulent Mach number, $\operatorname{\mathcal{M}}$, is a key ingredient for numerous star formation models (Krumholz & McKee, 2005; Federrath et al., 2010; Hennebelle et al., 2011; Federrath & Klessen, 2012; Federrath, 2013; Konstandin et al., 2016). We make the distinction between the scale-dependent turbulent Mach number,

[TABLE]

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,

[TABLE]

where $L$ is the cloud diameter, which corresponds to the outer scale of turbulence in our study (Federrath, 2013). For an isothermal cloud with purely turbulent dynamics, $\operatorname{\mathcal{M}}$ sets the width of the log-normal cloud density distribution,

[TABLE]

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 (Padoan et al., 1997; Passot & Vázquez-Semadeni, 1998; Kritsuk et al., 2007; Federrath et al., 2008, 2010; Konstandin et al., 2012b). The dispersion has been studied extensively and there have been many modifications to account for 2D projections of the 3D cloud (Burkhart & Lazarian, 2012), thermal and magnetic pressures (Padoan et al., 1997; Passot & Vázquez-Semadeni, 1998; Federrath et al., 2008; Price et al., 2011; Molina et al., 2012; Gazol & Kim, 2013), and non-isothermal (Nolan et al., 2015) and polytropic gases (Passot & Vázquez-Semadeni, 1998; Li et al., 2002; Federrath & Banerjee, 2015). 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 (Krumholz & McKee, 2005; Padoan & Nordlund, 2011; Hennebelle et al., 2011; Federrath & Klessen, 2012; Kainulainen et al., 2014). Beyond the density distribution, the turbulent Mach number may also play an important role in the distribution of cloud filament widths.

1.2 Filaments and the Turbulent Mach Number

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 (André et al., 2010; Men’shchikov et al., 2010; Schneider et al., 2012; André et al., 2014; Padoan et al., 2014; Arzoumanian, 2015; Federrath, 2016). 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 (Arzoumanian et al., 2011; Juvela et al., 2012; Palmeirim et al., 2013; André et al., 2014; Smith et al., 2014; Benedettini et al., 2015; Kirk et al., 2015; Federrath, 2016; Federrath et al., 2016; Smith et al., 2016; Arzoumanian et al., 2019). The standard deviation of the filament width distribution is thought to be associated with the sonic scale, $\ell_{s}$, in the clouds (Federrath et al., 2018). The sonic scale marks the transition between supersonic and subsonic velocity dispersions, and is theorised to be at length scale,

[TABLE]

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$ (Larson, 1981; Solomon et al., 1987; Ossenkopf & Mac Low, 2002; Heyer & Brunt, 2004; Kritsuk et al., 2007; Schmidt et al., 2009; Federrath et al., 2010; Roman-Duval et al., 2011; Federrath & Klessen, 2012; Federrath, 2016; Federrath et al., 2018). 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 (Federrath, 2016; Tokuda et al., 2018). However these are not the only structures that are formed and the densities of turbulent clouds have been shown to respect a fractal geometry.

1.3 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 (Falgarone & Phillips, 1996; Stutzki et al., 1998; Chappell & Scalo, 2001; Kauffmann et al., 2010; Schneider et al., 2013; Kainulainen et al., 2014; Rathborne et al., 2015). 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 (Scalo 1990; Elmegreen & Falgarone 1996; Sanchez et al. 2005; Kowal & Lazarian 2007; Federrath et al. 2009; Roman-Duval et al. 2010; Donovan Meyer et al. 2013; Konstandin et al. 2016; BFK19). However, where it falls between 2 and 3 depends upon the type of turbulent driving (Federrath et al., 2009), the rms $\operatorname{\mathcal{M}}$ (Konstandin et al. 2016; 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 $\operatorname{\mathcal{M}}$ limit) compared to 2D density slices ($\approx 1$ in the high $\operatorname{\mathcal{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 $\operatorname{\mathcal{M}}$, which is demonstrably an important quantity in star formation.

This study is organised into the following sections: In §2 we discuss the six cloud simulations that we use to construct our new Mach number - fractal dimension $(\operatorname{\mathscr{M}}-\operatorname{\mathcal{D}})$ relation. In §3 we summarise the fractal dimension method introduced by BFK19, including the key results. Next, in §4 we derive the new relation, and discuss the applications and limitations. Then in §5 we apply it to two quiescent subregions of the Polaris Flare to calculate the $\operatorname{\mathcal{M}}$ based purely on the fractal geometry of the column density. We compare this with previous estimates of the $\operatorname{\mathcal{M}}$ calculated in Schneider et al. (2013). Finally, in §6 we summarise our key findings.

2 Turbulent Molecular Cloud Models

In this study we use six purely hydrodynamical simulations of quiescent (non-star-forming) molecular clouds, with no self-gravitation, to construct our $\operatorname{\mathcal{M}}$ – $\operatorname{\mathcal{D}}$ relation. The parameter set of the molecular cloud models is listed in Table 1. For each of the simulations we solve the compressible Euler equations,

[TABLE]

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 Federrath et al. (2010), Federrath et al. (2018) and BFK19.

We run the six simulations with Mach numbers $\operatorname{\mathcal{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$ (Federrath et al., 2009; Price & Federrath, 2010). This gives us a wide set of $\operatorname{\mathcal{M}}$ values to construct the $\operatorname{\mathscr{M}}-\operatorname{\mathcal{D}}$ relation, encompassing transonic, slightly compressible flows, all the way to highly supersonic, and highly compressible flows that are saturated with shocks (Federrath, 2013). 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 1 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., Federrath et al. (2009); Federrath et al. (2010). This lets us perform our study only on the $xy$ projections (column densities), whilst still being representative 2D projections through any viewing angle.

3 Fractal Dimension Curves

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

3.1 Method Summary

There exists a power-law scaling between the mass and length scales (size) in real MCs (Larson, 1981; Myers, 1983; Falgarone & Phillips, 1996; Roman-Duval et al., 2010; Donovan Meyer et al., 2013). To calculate the mass-length dimension we utilise this power-law scaling, where the mass, $m$ is given by

[TABLE]

where we use the dimensionless $\ell/L$ for the length scales in the cloud, and $\operatorname{\mathcal{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^{\operatorname{\mathcal{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, $\operatorname{\mathcal{D}}$,

[TABLE]

We use Equation 8 to then define the fractal dimension at length scale $\ell_{i}/L$,

[TABLE]

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,\ldots,\ell_{i}/L\right\}$, where $\ell_{0}/L<\ell_{1}/L<\ldots<\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,

[TABLE]

to access the scaling in the column density, $\Sigma$.

To calculate the $\operatorname{\mathcal{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:

Identify the coordinates of the maximum column density pixel $\Sigma_{\text{max}}$, 2. 2.

expand a $(\ell/L)\times(\ell/L)$ square region centred on $\Sigma_{\text{max}}$, creating our length scale hierarchy, 3. 3.

calculate the mass $m$, within each of the $(\ell/L)\times(\ell/L)$ squares, 4. 4.

use the relation shown in Equation 9 to determine the fractal dimension, $\operatorname{\mathcal{D}}$, as a function of $\ell_{i}/L$, always using all length scales below $\ell_{i}/L$ to calculate $\operatorname{\mathcal{D}}$ on $\ell_{i}/L$.

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.

3.2 Key Results from the $\operatorname{\mathcal{D}}(\ell)$ Curves

In BFK19 we show that the fractal dimension curves from each simulation (shown in Figure 2) 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 $\operatorname{\mathcal{D}}(\ell/L)$, which models a smooth transition between space-filling clouds ($\operatorname{\mathcal{D}}=2$ in the 2D projection) and clouds saturated with planar shocks ($\operatorname{\mathcal{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,

[TABLE]

where $\operatorname{erfc}$ is the complimentary error function, $\operatorname{\mathcal{D}_{\text{min}}}$ and $\operatorname{\mathcal{D}_{\text{max}}}$ 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 2. This fit encodes the limits $\lim_{(\ell/L)\rightarrow 1}\operatorname{\mathcal{D}}(\ell/L)=\operatorname{\mathcal{D}_{\text{min}}}$ and $\lim_{(\ell/L)\rightarrow 0}\operatorname{\mathcal{D}}(\ell/L)=\operatorname{\mathcal{D}_{\text{max}}}$ 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 $\operatorname{\mathcal{D}}(\ell/L)$ are shown in the sub-panel of Figure 3. The plot shows that the $\operatorname{\mathcal{D}}$ of the column density is bounded between $\operatorname{\mathcal{D}_{\text{max}}}=2$ and $\operatorname{\mathcal{D}_{\text{min}}}=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 $\operatorname{\mathscr{M}}-\operatorname{\mathcal{D}}$ relation.

4 The $\operatorname{\mathscr{M}}$ – $\operatorname{\mathcal{D}}$ Relation

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 1. This link is made available to us by approximating $\operatorname{\mathscr{M}}(\ell)$ using scaling relations from models of supersonic turbulence.

4.1 Constructing the Relation

Using the second-order structure function, $\operatorname{\text{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\ldots\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 Konstandin et al. 2012a). This construction of the Mach number follows a power-law of the form,

[TABLE]

where $p\approx 1/2$ for supersonic turbulence and $p\approx 1/3$ for subsonic turbulence (Kolmogorov, 1941; Burgers, 1948; Kritsuk et al., 2007; Schmidt et al., 2009; Konstandin et al., 2012a; Federrath, 2013; Federrath et al., 2018). We use the $p=1/2$ case, $\ell/L\sim\operatorname{\mathscr{M}}(\ell/L)^{2}$, to transform all length scales into Mach numbers. We set $\ell=L$ to find the constant of proportionality, $\operatorname{\mathcal{M}}^{-2}$, i.e. on large scales in the cloud $\operatorname{\mathscr{M}}$ is $\approx\operatorname{\mathcal{M}}$ (Federrath, 2013). Hence the transformation is

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

and

[TABLE]

This corresponds to the turbulent Mach number on the length scale that $\operatorname{\mathcal{D}}$ was measured, since $\operatorname{\mathscr{M}}(\operatorname{\mathcal{D}}(\ell))$. The values of the estimated and derived parameters for this fit are shown in Table 3.

4.2 $\operatorname{\mathscr{M}}$ – $\operatorname{\mathcal{D}}$ Relation Results

In the main plot of Figure 3 we show $\operatorname{\mathscr{M}}$ as a function of $\operatorname{\mathcal{D}}$, derived from the column density data. The black line shows Equation 16, 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 $\operatorname{\mathcal{M}}$ the relation is extremely sensitive to small changes in $\operatorname{\mathcal{D}}$. The second is that for high $\operatorname{\mathcal{M}}$ the temporal fluctuations of $\operatorname{\mathcal{D}}$ become significant, spanning over $>0.2\operatorname{\mathcal{D}}$. This means that our relation will perform best at measuring regions of MCs with $10\gtrsim\operatorname{\mathcal{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 3 may only work well, without modification, on quiescent clouds, without significant deviation from natural mixing, since $\operatorname{\mathcal{D}}$ is sensitive to changes in driving Federrath et al. (2009). Acknowledging these limitations, we now test the performance of the relation on $Herschel$ observations of subregions from the Polaris Flare cloud.

5 Application on quiescent subregions in the Polaris Flare Cloud

The Polaris Flare is a high Galactic latitude cloud which is located at a distance $\leq 150$pc (Falgarone et al., 1998; Miville-Deschênes et al., 2010; Schlafly et al., 2014). It has weak, but significant CO emissions in regions of higher hydrogen column density (Falgarone et al., 1998; Meyerdierks & Heithausen, 1996; Miville-Deschênes et al., 2010). There is no active star-formation in Polaris, only 5 starless cores were detected (André et al., 2010; Ward-Thompson et al., 2010) 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.

5.1 Observational Data

The Polaris region was observed as part of the $Herschel$ Gould Belt survey (HGBS, André et al. 2010) using the PACS and SPIRE instruments on-board $Herschel$. For all observational details, we refer to (André et al., 2010; Men’shchikov et al., 2010; Ward-Thompson et al., 2010; Miville-Deschênes et al., 2010). 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 Bernard et al. (2010). Column density and temperature maps were then produced at an angular resolution of 18*′′, following the procedure outlined in Palmeirim et al. (2013) 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}$ cm2 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 André et al. (2010) and Könyves et al. (2015) 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 (Miville-Deschênes et al., 2010; Schneider et al., 2013). 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 $\operatorname{\mathcal{M}}$ in Polaris (Schneider et al., 2013), the ‘saxophone’ and ‘quiet’ subregions, shown in Figure 4. Schneider et al. (2013) made estimates of the turbulent Mach number for the two regions based on the CO111The CO data stem from 12CO 2$\to$1 and 13CO 1$\to$0 observations from the KOSMA and FCRAO telescopes published in Bensch et al. (2003). velocity dispersion. The authors calculated the $\operatorname{\mathcal{M}}$ for each subregion using,

[TABLE]

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 $\operatorname{\mathcal{M}}$ are shown in column (2) of Table 4.

5.2 Application and considerations

Using the method summarised in §3.1 we construct the fractal dimension curves for each of the two regions and then use Equation 16 to calculate the Mach number222Our fractal dimension implementation on the Polaris Flare cloud is available here: 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 $\operatorname{\mathcal{D}}$ when $\operatorname{\mathcal{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 $\operatorname{\mathcal{D}}$, as described in §3.1. We do this because it allows us to calculate $\operatorname{\mathcal{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 §1.1 and 1.2, respectively). It also allows for comparison with previous estimates made for $\operatorname{\mathcal{M}}$ in the ‘quiet’ and ‘saxophone’ subregions of the Polaris Flare molecular cloud using Equation 17.

5.3 $\operatorname{\mathcal{M}}(\operatorname{\mathcal{D}})$ Results

In Figure 5 we show our $\operatorname{\mathcal{M}}$ estimates for the two subregions of Polaris, indicated in red, and compare them with the $\operatorname{\mathcal{M}}$ estimated by Schneider et al. (2013), shown as blue-dashed regions. The $\operatorname{\mathcal{M}}$ calculated by Schneider et al. (2013) and the $\operatorname{\mathcal{M}}$ calculated in this study are shown in Table 4. We find $\operatorname{\mathcal{M}}=2\pm 1$ and $\operatorname{\mathcal{M}}=10^{+55}_{-4}$, whereas Schneider et al. (2013) finds $\operatorname{\mathcal{M}}3\pm 1$ and $\operatorname{\mathcal{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 $\operatorname{\mathcal{D}}$, especially for $\operatorname{\mathcal{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 4. To understand why there may be differences between the $\operatorname{\mathcal{M}}(\operatorname{\mathcal{D}})$ and previous estimates based on the CO velocity dispersion we turn to the calculation of $\operatorname{\mathcal{D}}$.

To estimate $\operatorname{\mathcal{M}}$ we first calculate $\operatorname{\mathcal{D}}$, which is shown in column (3) of Table 4. We find ‘quiet’, the lower $\operatorname{\mathcal{M}}$ region, has a $\operatorname{\mathcal{D}}=1.76\pm 0.05$, and ‘saxophone’, the higher $\operatorname{\mathcal{M}}$ region, $\operatorname{\mathcal{D}}=1.60\pm 0.04$. This is consistent with BFK19, who argues that for higher $\operatorname{\mathcal{M}}$ flows we should expect lower $\operatorname{\mathcal{D}}$, corresponding to the introduction of compressive shocks into a diffuse cloud with increasing $\operatorname{\mathcal{M}}$, and previous studies have calculated $\operatorname{\mathcal{D}}=1.4-1.8$ for column densities (Elmegreen & Falgarone, 1996; Elmegreen & Scalo, 2004; Sanchez et al., 2005; Rathborne et al., 2015). This suggests that the $\operatorname{\mathcal{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 4), which has a density-length scaling relation $\rho\sim\ell^{-2}$, or $\Sigma\sim\ell^{-1}$, and will act to reduce $\operatorname{\mathcal{D}}$ in the vicinity of $\Sigma_{\text{max}}$ (Schneider et al., 2013; Federrath, 2016; André, 2017). This may account for the slightly higher $\operatorname{\mathcal{M}}$ that we estimate for ‘saxophone’, compared to the Mach estimate in Schneider et al. (2013). For the ‘quiet’ subregion we slightly underestimate $\operatorname{\mathcal{M}}$. This is because the column density (see the right column density map in Figure 4) is diffuse, and lacks the shock structures that we see introduced between the $\operatorname{\mathcal{M}}=1$ and $\operatorname{\mathcal{M}}=4$ simulations in Figure 1. 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\operatorname{\mathcal{D}}$ for the mass-length method (Federrath et al., 2009).

6 Summary and Key Findings

In this study we construct a new empirical relation for the scale-dependent three-dimensional (3D) Mach number, $\operatorname{\mathscr{M}}$, and the fractal dimension, $\operatorname{\mathcal{D}}$, of the column density for turbulent clouds. We use the mass-length fractal dimension method introduced in Beattie et al. (2019) (BFK19) on six hydrodynamical cloud simulations, with root-mean-squared (rms) Mach number, $\operatorname{\mathcal{M}}$, varying from $\operatorname{\mathcal{M}}=1$ to $100$. We apply the method on the column densities, with an example of the densities shown Figure 1, to construct $\operatorname{\mathcal{D}}$ as a function of length scale, $\ell$. We then transform the cloud length scales to $\operatorname{\mathscr{M}}$ using the scaling relation $\ell\sim\operatorname{\mathscr{M}}^{2}$, for supersonic turbulence (Burgers, 1948; Federrath, 2013). Using this data we are able to construct $\operatorname{\mathcal{D}}(\operatorname{\mathscr{M}})$, and finally $\operatorname{\mathscr{M}}(\operatorname{\mathcal{D}})$. Using $\operatorname{\mathscr{M}}(\operatorname{\mathcal{D}})$ we construct $\operatorname{\mathcal{M}}$, where $\operatorname{\mathcal{M}}\sim\operatorname{\mathscr{M}}(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 Schneider et al. (2013). We summarise our key findings below:

• •

We propose a new empirical relation for the scale-dependent Mach number and the fractal dimension of the column density,

[TABLE]

as defined in Equation 16 and plotted in Figure 3, where $\xi_{1}=0.9\pm 0.1$, $\xi_{2}=0.2\pm 0.2$, $\Omega=0.22\pm 0.07$ and $\operatorname{\mathcal{D}_{\text{min}}}$, the fractal dimension of the column density in the high $\operatorname{\mathcal{M}}$ limit is $1.55\pm 0.13$. This relation allows for the calculation of $\operatorname{\mathcal{M}}$ for clouds in the range $10\gtrsim\operatorname{\mathcal{M}}\gtrsim 0.1$. Very large and very low $\operatorname{\mathcal{M}}$ are inappropriate for the model due to the steep tails in the inverse complimentary error function.

• •

We use the mass-length fractal dimension method in BFK19 to calculate the $\operatorname{\mathcal{D}}$ of the ‘saxophone’ and ‘quiet’ subregions of the Polaris Flare, shown in Figure 4. We find $\operatorname{\mathcal{D}}=1.60\pm 0.04$ and $\operatorname{\mathcal{D}}=1.76\pm 0.05$ for ‘saxophone’ and ‘quiet’, respectively, consistent with the thesis that higher Mach number flows reduce $\operatorname{\mathcal{D}}$, by turning diffuse, space-filling structures into compressive shocks and filaments.

• •

Using $\operatorname{\mathcal{D}}$ we estimate $\operatorname{\mathcal{M}}$ for each of the subregions. We find ‘quiet’ has a $\operatorname{\mathcal{M}}\sim 2$ and ‘saxophone’ has a $\operatorname{\mathcal{M}}\sim 10$, shown in Figure 5. This is comparable to the estimates made in Schneider et al. (2013), $\operatorname{\mathcal{M}}\sim 3$ and $\operatorname{\mathcal{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 16, 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 $\operatorname{\mathcal{D}}$.

• •

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.

Acknowledgements

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).