© ESO, 2017

1. Introduction

In the prevalent paradigm of the turbulence-regulated interstellar medium (ISM) and star formation, the star formation rate of the ISM is linked to the internal density distribution of molecular clouds. In this picture, the cloud structure is strongly affected by supersonic turbulent motions that drive the formation of density enhancements, some of which become self-gravitating and form stars (for reviews, see Elmegreen & Scalo 2004; Hennebelle & Falgarone 2012; Padoan et al. 2014). Specifically, theoretical works predict, and star formation rate models assume, that the initial density distribution of the ISM is characterised by a log-normal probability density function of gas densities (hereafter ρ-PDF), a form that naturally results in isothermal, non-gravitating, supersonically turbulent gas (e.g. Vazquez-Semadeni 1994; Padoan et al. 1997; Passot & Vázquez-Semadeni 1998). The width of the lognormal ρ-PDF, σ_s, where s = lnρ/ρ₀ and ρ₀ is the mean density, depends on turbulence energy, magnetic field strength, and the fraction of compressive energy in the gas (Molina et al. 2012; Federrath & Banerjee 2015) ${\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\ln (\mathrm{1}\mathrm{+}{\mathit{b}}^{\mathrm{2}}{\mathrm{ℳ}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathit{\beta }}{\mathrm{1}\mathrm{+}\mathit{\beta }})\mathit{,}$(1)where ℳ_s is the sonic Mach number, β is the ratio between thermal and magnetic energies, and b describes the mixture of solenoidal (divergence-free) and compressive (curl-free) modes in the acceleration field that drives the turbulence (Federrath et al. 2008, 2010). Equation (1) represents a fundamental prediction of the turbulence-regulated ISM paradigm, and it plays a crucial role in analytic theories of star formation by coupling the density statistics of the gas to the physical processes acting in it (e.g. Krumholz & McKee 2005; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011; Federrath & Klessen 2012, for a review, see Padoan et al. 2014).

However, few works have obtained observational constraints for the relationship given by Eq. (1). This is because of the observational difficulties in probing it: none of its parameters are directly accessible to observations. The sonic Mach number can be estimated through measurements of the molecular line widths, but it is not clear what fraction of the observed line width is caused by turbulence. The ρ-PDF width cannot be directly measured because observations only probe the projected column densities, not volume densities. In addition, observational techniques all probe only a limited range of column densities (e.g. Goodman et al. 2009), making a reliable quantification difficult. It is also difficult to measure magnetic field strengths in clouds, leaving the β parameter uncertain. Consequently, the relationship given by Eq. (1) remains poorly constrained.

Works targeting single molecular clouds in the solar neighbourhood have suggested values of b ~ 0.5 (Padoan et al. 1997; Brunt et al. 2010). Ginsburg et al. (2013) found that turbulence in the massive star formation region W49 is likely driven compressively with b > 0.4. Using small samples of clouds to probe the mean value of b, Kainulainen & Tan (2013) derived $\mathit{b}\mathrm{=}\mathrm{0.}{\mathrm{2}}_{-0.22}^{\mathrm{+}\mathrm{0.37}}$, assuming non-magnetised gas. Kainulainen et al. (2013) compared the dense gas mass fractions of molecular clouds with simulations that spanned a variety of parameters {b,β, ℳ_s} and concluded that only simulations with b ≈ 1/3−0.4 match the observations of the solar neighbourhood clouds. Kainulainen et al. (2013) also found that infrared dark clouds possibly showed higher b values. The above works tend to suggest dominance of compressive driving modes (b ≳ 0.4) in the clouds located in the Galactic disc and spiral arms. However, Federrath et al. (2016) found b = 0.22 ± 0.12 in the central molecular zone cloud G0.253+0.016 (“the Brick”), which indicates that the driving is primarily solenoidal in environments that are dominated by strong shearing motions. In summary, only a few observational works have constrained b so far, and they are based on individual measurements or small samples and employ different approaches to estimate b; systematic studies are lacking to date.

Meanwhile, numerical simulations of molecular cloud formation in disc galaxies and idealised colliding flow setups predict a range of b values emerging in different cloud environments and at different times (Jin et al. 2017; Körtgen et al. 2017). These studies suggest that a single constant b value cannot be used to describe the turbulence driving of all clouds, but instead that b can vary significantly from cloud to cloud, covering the full range from solenoidal to compressive driving.

Enabling progress in the topic, we have developed a technique to determine the volumetric PDF widths (σ_s) from projected (column) density maps. We have previously used the technique to perform the first systematic observational quantification of the PDF widths in the solar neighbourhood clouds (Kainulainen et al. 2014). In this paper, we analyse the ℳ_s – σ_s relationship in the solar neighbourhood clouds using the σ_s data from Kainulainen et al. (2014) and ℳ_s data from the literature. We also derive new σ_s values for the Orion A, Orion B, and California molecular clouds. This allows us to compare Eq. (1) with the latest observational results and to provide constraints for the analytic star formation rate models using Eq. (1)

2. Data

We adopted most data for σ_s and ℳ_s from the literature. Kainulainen et al. (2014) derived σ_s values for several solar neighbourhood clouds, and we use these values (listed in Table 1). In short, their technique was based on reconstructing the observed column density data with the help of an ensemble of three-dimensional model forms (prolate spheroids). The model forms are arranged hierarchically, nested inside each other, which allows construction of complex structures, such as fractals, which contain filaments and sheets, similar to observed cloud structures. The technique was used in conjunction with near-infrared dust extinction derived column density maps from Kainulainen et al. (2009) that had a spatial resolution of 0.1 pc, and they probed the column densities between N(H₂) ≈ 1−25 × 10²¹ cm^-2. Kainulainen et al. (2014) estimated the uncertainty of the resulting σ_s values to be about 20%. We refer to Kainulainen et al. (2014) for further details.

We adopted ℳ_s values for most clouds from Kainulainen & Tan (2013), who used CO (1-0) line emission data from Dame et al. (2001) to estimate the total line widths of the clouds (listed in Table 1). The sonic Mach numbers were computed from the line widths using ${\mathrm{ℳ}}_{\mathrm{s}}\mathrm{=}\frac{\sqrt{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{1}\mathrm{D}}}{{\mathit{c}}_{\mathrm{s}}}\mathit{,}$(2)where ${\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{1}\mathrm{D}}$ is the observed line width and c_s the isothermal sound speed at 15 K. The Musca cloud is not covered by Dame et al. (2001), and we adopted the Mach number for it from Hacar et al. (2016). We adopted 5 K as the 3σ uncertainty in the temperature, resulting in a 22% uncertainty (3σ) in the sonic Mach number. In addition to the temperature uncertainty, it is unclear how well the Mach number derived from CO reflects the turbulent energy of the cloud. Szűcs et al. (2016) analysed numerical simulations that included gas-phase chemistry and radiative heating and cooling, and found that the observed CO line widths are within 30–40% of the true velocity dispersion. Taking this into account, we assigned a total uncertainty of 60% (3σ) to the Mach number.

To expand the sample provided by the above works, we additionally derived ρ-PDFs for Orion A, Orion B, and the California molecular clouds using the technique of Kainulainen et al. (2014). We summarise the derivation and present the resulting ρ-PDFs in Appendix A. We adopted the CO line widths for Orion A, Orion B, and the California cloud from Dame et al. (2001).

In summary, our sample contains 15 solar neighbourhood molecular clouds and includes most of the major cloud complexes within the distance of about 450 pc. The relevant parameters of the clouds are listed in Table 1.

3. Results

Figure 1 presents the relationship between the sonic Mach numbers and density distribution widths in the 15 solar neighbourhood clouds. The ℳ_s values span the range [4, 17] and the σ_s values the range [1.29, 2.08]. It is not trivial to assess whether the variables are correlated because the probability density functions of the uncertainties in the variables are not well known. The possibility of having a correlation can be probed by various tests. For example, ignoring all uncertainties, the Pearson correlation coefficient for the data is 0.35 and the p-value is 0.20, which indicates no correlation. Similarly, a linear fit using errors in both variables (the fitexy routine in IDL; Press & Teukolsky 1992) results in a slope of 0.02 ± 0.02, indicating no significant correlation. However, we cannot rule out a weak correlation either (discussed further below). Clearly, including the sample clouds over a wider range of Mach numbers, especially including extreme ℳ_s values, would be beneficial in further studies; we show in Appendix B that a sample of 50–100 clouds is sufficient to reliably establish a correlation, if present.

We next considered the observed ℳ_s – σ_s data in the context of turbulence-regulated density variance. Figure 1 shows the prediction given by Eq. (1) for the hydrodynamic case (β = ∞) and for a moderately magnetised case (β = 0.3, that is, B ≈ 4μG at n = 100 cm^-3 and T = 15 K). The comparison of these models with the observed data gives rise to two main results.

First, the data points agree well in general with the ranges set by the models. The spread of the models that have no magnetic field (b varies) covers all but one data point. If the magnetic field is stronger than given by β ≲ 0.05 (B ≳ 30μG at n = 1000 cm^-3 and T = 15 K), all but two clouds have σ_s values in disagreement with the models. This implies a general agreement between the turbulence-regulated density structure and observations if the magnetic field is relatively weak.

Second, given the large uncertainties, the data cannot rule out a weak correlation such as is predicted by Eq. (1); the correlation coefficient of the observed data is in agreement with a sample drawn from the model represented by Eq. (1) (demonstrated further in Appendix B). Furthermore, the large uncertainties enable explaining all the observations with a single ℳ_s – σ_s relationship (see Appendix B). A single ℳ_s– σ_s relationship corresponds to a family of parameters {b,β} (see Eq. (1)). This family can be characterised by a parameter a = b²β/ (1 + β). We determined the best-fitting a parameter by fitting Eq. (1) to the observed data with a as the free parameter. This yielded the best-fitting model with a = 0.30 ± 0.06 (shown in Fig. 1, top panel, as the dotted line). The constraints on the above results present for b and β are discussed in the following Sect. 4.

4. Discussion

The interpretation of the observed scatter in the ℳ_s – σ_s data remains uncertain because the data are limited. All clouds can share one set of { β,b } values, or they can have individually varying β and/or b values; the observations are in agreement with either case. Distinguishing between these two options would require more accurate measurements of ℳ_s and σ_s  or alternatively, a significantly larger sample of clouds. We next discuss how b and β can be constrained in the above two cases.

In the former case, that is, if all clouds originated from one { β,b } family, we can use the best-fitting model to constrain β and b (Fig. 1). The best-fitting model allows us to rule out several regimes: a regime of strong compression combined with a low magnetic field (b ≳ 0.7 and β ≳ 3), weak compression (b ≲ 0.35), or a strong magnetic field (β ≲ 0.1).

If the { β,b } values vary from cloud to cloud, we can obtain estimates for b of each cloud by fixing β. Lower limits for b result from the hydrodynamic case, that is, β = ∞ (Table 1). The lower limits are in agreement with the picture of turbulence-regulated density variance for 13 out of 15 clouds, which means that they are between 0.3 and 1. The median lower limit for b is 0.54.

Estimates of b beyond the lower limits can be achieved when the magnetic field strengths are known. However, such measurements are very difficult to obtain. The most direct measurement is provided by the Zeeman-splitting measurements, which indicate a maximum line-of-sight field strength of about 10 μG at densities lower than 1000 cm^-3 (Crutcher 2012). Alternatively, the mean magnetic field strengths for several nearby clouds have been estimated from the Planck data using the Davis-Chandrasekhar-Fermi technique (CDF, Planck Collaboration Int. XXXV 2016). The technique results in upper limits for the field strengths (see the discussion in Planck Collaboration Int. XXXV 2016). The field strengths reported in Planck Collaboration Int. XXXV (2016) are 10–50 μG, in reasonable agreement with, but also systematically higher than, the results from Zeeman-splitting. These field strengths correspond to a lower-limit range for β ~ 0.002–0.05 in the density of 100 cm^-3 and to a temperature of 15 K.

Given the uncertainty of the magnetic field strength estimates, we performed a simple exercise. We computed for each cloud a b_MHD value assuming a field strength of 8 μG (see Table 1), broadly in agreement with Zeeman-splitting measurements and the upper limits from CDF estimates (Crutcher 2012; Planck Collaboration Int. XXXV 2016). At the ρ = 100 cm^-3 and T = 15 K, β = 0.08 follows. Almost all b_MHD values exceed unity, which is in contradiction with the framework of Eq. (1) (values between ~0.3–1 are expected, Federrath et al. 2008, 2010). Low b values are not in agreement with observations; 3 clouds out of 15 have b(β = ∞) below 0.4 and none have b_MHD below 0.4.

Several issues can contribute to the contradiction between the b_MHD values and Eq. (1). Federrath et al. (2016) argued that the turbulent magnetic field component, not the total field strength, is relevant for Eq. (1). This is because in the derivation of Eq. (1), the magnetic pressure in the post-shock region enters (Molina et al. 2012), which is likely dominated by the turbulent field component instead of the total field strength. The turbulent field component is often significantly smaller than the total field strength, especially in the presence of strong guide fields (Pillai et al. 2015; Federrath et al. 2016; Federrath 2016). Using the turbulent field component (instead of the total field component) would decrease the b_MHD values, but the turbulent field component is difficult to measure observationally. Observational estimates of the magnetic field strength are indeed usually very uncertain; it may be that our adopted values of 8 μG at 100 cm^-3 are not good estimates. Choosing 5 μG at 200 cm^-3 instead leads to β = 0.41 and to majority of the b_mhd values being lower than or equal to unity (values given in Table 1). It may also be that the density variance suffers from systematics; the employed technique was tested against numerical simulations (see Kainulainen et al. 2014), but those simulations do not necessarily capture all relevant aspects of real molecular clouds. Finally, the correspondence between ℳ_s and turbulence energy is also unclear and could be tested in more detail with numerical simulations (e.g. Szűcs et al. 2016).

Fig. 1 Top: relationship between the density distribution width, σs, and sonic Mach number, ℳs. The red dashed curves show the case β = ∞, the blue curves β = 0.3, and green dash-dotted curves β = 0.05. The thicker (and higher) curves correspond to b = 1 and the thinner (and lower) curves to b = 1/3. The dotted line shows the best-fitting model with b2β/ (1 + β) = 0.30 ± 0.06. Bottom: the relationship between the magnetic field strength, β, and the compression parameter, b, for the model that best fits all data points. The thin solid, dashed, and dotted lines show the 1, 2, and 3σ confidence levels.

5. Conclusions

We analysed the relationship between the ρ-PDF width, σ_s, and sonic Mach number, ℳ_s, in 15 solar neighbourhood molecular clouds. We used data primarily from the literature, but also derived new ρ-PDFs for the Orion A, Orion B, and California molecular clouds. Our conclusions are as follows:

• 1.

  The 15 clouds span the ranges ℳ_s = [4, 17] and σ_s = [1.29, 2.08]. The variables show no significant correlation, but because of the large scatter in the data, we cannot rule out a weak correlation either.

• 2.

  In the context of turbulence-dominated gas, the observed ℳ_s – σ_s relationship can constrain the turbulence-driving parameter b and the magnetic field strength measured by the thermal-to-magnetic pressure ratio, β. The majority of the clouds is not in agreement with field strengths stronger than β < 0.05, and none with β < 0.01. When we consider the model that best fits the 15 clouds together, the data rule out the regime of strong compression combined with a weak magnetic field: b ≳ 0.7 and β ≳ 3. The data also rule out weak compression, b ≲ 0.35, or a strong magnetic field, β ≲ 0.1.

• 3.

  When combined with magnetic field strength estimates, the ℳ_s – σ_s data rule out solenoidal driving with b < 0.4. However, most clouds have compression parameters larger than unity, which is not expected in the context of turbulence-regulated density variance (Eq. (1)). It is not trivial to asses the source of the contradiction; one likely explanation is that it is not the total magnetic field strength, but the smaller-scale turbulent magnetic field component (i.e. the turbulent magnetic pressure in the post-shock regions) that affects the density variance based on the derivation of Eq. (1) by Padoan & Nordlund (2011) and Molina et al. (2012), as discussed in Federrath et al. (2016).

Observational uncertainties restrict the determination of the parameters entering Eq. (1), but we demonstrated that systematic studies can start placing constraints on them; using our techniques, a sample of 50–100 clouds can establish the correlation between σ_s and ℳ_s, if present (Appendix B).

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 639459 (PROMISE). C.F. acknowledges funding provided by the Australian Research Council’s Discovery Projects (grants DP150104329 and DP170100603) and by the ANU Futures Scheme.

References

Appendix A: ρ-PDFs of Orion A, Orion B, and the California cloud

We derive ρ-PDFs for the Orion A, Orion B, and the California molecular clouds using the three-dimensional modelling technique presented by Kainulainen et al. (2014, see the paper for the details of the technique). As the data for the technique, we used for Orion A and the California nebula Herschel dust emission derived column density data from Stutz & Kainulainen (2015) and Harvey et al. (2013), respectively. These data have a spatial resolution of a full-width at half-maximum (FWHM) = 38′′, which corresponds to about 0.08 pc at the distance of Orion and the California nebula. This equals the physical resolution of the data used to derive the other ρ-PDFs we employed in this work (see Kainulainen et al. 2009). For Orion B, we used a near-infrared extinction map from Kainulainen et al. (2009) that has a resolution of FWHM = 150′′. Acknowledging the caveat that this is coarser than the data for other clouds, we include the cloud in the sample. We tested the possible effect of the resolution by smoothing the Orion A data with a factor of three and repeating the ρ-PDF. This did not change the resulting ρ-PDF width, although it did truncate the ρ-PDF as the highest column densities were smoothed out. The final ρ-PDFs or Orion A, Orion B, and the California nebula are shown in Fig. A.1 and the resulting ρ-PDF widths are listed in Table 1.

Fig. A.1 ρ-PDFs derived for Orion A (top), Orion B (middle), and the California molecular clouds (bottom). The solid line is a fit of a log-normal function to the data.

Appendix B: On the scatter of the observed ℳ_s – σ_s relationship

We analyse here the scatter in the observed ℳ_s – σ_s relationship to address three questions: 1) can the data be well fitted by a single model represented by Eq. (1); 2) can the observations rule out a correlation between ℳ_s and σ_s such as is predicted by Eq. (1); and 3) given the uncertainties, can a larger sample reliably establish the presence of a relationship described by Eq. (1), if present?

We address the first two questions with a Monte Carlo simulation in which we drew random samples of 15 data points from the one model (Eq. (1)) that best fits all the observed data points (see Fig. 1). The model is characterized by the parameter a = 0.30 ± 0.06. We restricted the draw between ℳ_s = 4–20. We drew 15 {ℳ_s, σ_s} pairs from the model so that the probability P along the model curve, that is, dP/ df where $\mathit{f}\mathrm{=}\mathrm{(}\ln \mathrm{1}\mathrm{+}{\mathit{b}}^{\mathrm{2}}{\mathrm{ℳ}}_{\mathrm{s}}^{\mathrm{2}}\mathit{\beta }\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, is constant and uniform. The data points were then assigned an error, according to the observational uncertainties, that is, 20% for both σ_s and ℳ_s. This procedure was repeated 10⁴ times to obtain reasonable statistics.

For each repetition, we computed the reduced chi-square value between the model and the simulated data points. The chi-square distribution is shown in Fig. B.1 together with the chi-square of the observations and the best-fitting model. It shows that the observed chi-square does not strongly deviate from the distribution of the chi-square of the simulation. This indicates that the observed ℳ_s – σ_s relationship can be explained with a single model represented by Eq. (1).

For each repetition, we also computed the Pearson correlation coefficient. The resulting distribution of coefficients is shown in Fig. B.1. For reference, a correlation coefficient higher than ~0.7 is usually taken to indicate a significant correlation (p < 0.003). The observed correlation coefficient is well within the range of the simulated coefficients. Importantly, the distribution shows that in ~8% of the cases, a significant correlation is detected. This indicates that we cannot rule out a correlation such as predicted by Eq. (1) with our data, and in fact, with only 15 data points, we would not expect to be able to do so (given the uncertainties of the measurements).

Finally, we address the third question with the following exemplary setup. An observational sample that has the best chance of detecting a correlation between σ_s and ℳ_s, if present, is one that probes extreme Mach numbers, both small and large. With that in mind, we considered a sample of N_clouds, half of which were chosen to lie between ℳ_s= [5,10] and the other half between [15,20]. How large is the sample that is required to detect a correlation, if any, assuming the uncertainties of our current observational sample? To determine this, we repeated the Monte Carlo simulation described above, but changed the number of data points, N_clouds, between 15 and 200, and restricted the Mach numbers at the ranges 5–10 and 15–20. For each simulation characterised by N_clouds, we saved the correlation coefficient and computed the corresponding p-value. Figure B.1 shows as an example histograms for N_clouds = 15,30,50, and 80. The experiment shows that as the sample size grows, the ability to detect the correlation of the variables improves. With 80 clouds, the fraction of p-values above 0.003 is smaller than 0.5%. Even though this experiment is simplistic, it indicates that using our adopted techniques, a well-constructed sample of some 50–100 clouds will be sufficient to detect the correlation between σ_s and ℳ_s, if present.

Fig. B.1 Top: distribution of the reduced chi-square values in the simulation. The solid vertical lines show the values from the observed data. Middle: the distribution of the Pearson correlation coefficients in the simulation. The limit for significant correlation is about 0.7. The solid vertical lines show the values from the observed data. Bottom: distributions of p-values in a simulation with 15 (red), 30 (purple), 50 (blue), and 80 (green) clouds. The dotted line shows the p-value 0.003.

All Tables

All Figures

Fig. 1 Top: relationship between the density distribution width, σs, and sonic Mach number, ℳs. The red dashed curves show the case β = ∞, the blue curves β = 0.3, and green dash-dotted curves β = 0.05. The thicker (and higher) curves correspond to b = 1 and the thinner (and lower) curves to b = 1/3. The dotted line shows the best-fitting model with b2β/ (1 + β) = 0.30 ± 0.06. Bottom: the relationship between the magnetic field strength, β, and the compression parameter, b, for the model that best fits all data points. The thin solid, dashed, and dotted lines show the 1, 2, and 3σ confidence levels.

In the text

	Fig. A.1 ρ-PDFs derived for Orion A (top), Orion B (middle), and the California molecular clouds (bottom). The solid line is a fit of a log-normal function to the data.
In the text

Fig. B.1 Top: distribution of the reduced chi-square values in the simulation. The solid vertical lines show the values from the observed data. Middle: the distribution of the Pearson correlation coefficients in the simulation. The limit for significant correlation is about 0.7. The solid vertical lines show the values from the observed data. Bottom: distributions of p-values in a simulation with 15 (red), 30 (purple), 50 (blue), and 80 (green) clouds. The dotted line shows the p-value 0.003.

In the text