The Transition from a Lognormal to a Power-law Column Density Distribution in Molecular Clouds: An Imprint of the Initial Magnetic Field and Turbulence

The column density probability distribution function (PDF) provides an effective way to analyze the dynamics and the evolution of molecular clouds from both observational (e.g., Burkhart et al. 2015b; Schneider et al. 2015a, 2016; Pokhrel et al. 2016) and theoretical (e.g., Burkhart 2018; Körtgen et al. 2019) perspectives. Numerical simulations have established that non-self-gravitating gas with driven turbulence results in many interacting shocks that yield a lognormal density or column density PDF (e.g., Padoan et al. 1997; Passot & Vázquez-Semadeni 1998; Scalo et al. 1998; Federrath et al. 2008; Molina et al. 2012). The addition of self-gravity into simulations introduces a high-density power-law tail to the PDF (Ballesteros-Paredes et al. 2011; Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2013; Ward et al. 2014; Auddy et al. 2018). Observed column density PDFs have an underlying lognormal shape with an additional power-law tail (Kainulainen et al. 2009; Alves et al. 2014) that starts at a transitional column density (Schneider et al. 2015b).

A lognormal shape is associated with quiescent clouds that do not have active star formation (Kainulainen et al. 2009; Lombardi et al. 2015; Schneider et al. 2015b). In contrast, the active star-forming clouds have an excess of high column density with a prominent power-law tail as well as a lognormal peak. The lognormal feature is often considered to be a direct imprint of driven supersonic turbulence and its width is attributed to the strength of the sonic Mach number (Collins et al. 2012; Molina et al. 2012; Burkhart et al. 2015a). However, it has also been suggested that the lognormal may be a more general characteristic that is set by both supersonic turbulence and gravitationally driven ambipolar diffusion (Tassis et al. 2010) or global gravitational contraction (Ballesteros-Paredes et al. 2011).

The power-law part of the PDF (Schneider et al. 2013; Alves et al. 2017) is a signature of gravitational contraction (due to the self-gravity of the gas) and can be associated with the formation of condensed cores. The PDF ${dN}/d\mathrm{log}{\rm{\Sigma }}\propto {{\rm{\Sigma }}}^{-\alpha }$ has an index α = 2 in the limit of isothermal gravitational contraction (see, e.g., Appendix A of Auddy et al. 2018), and is set by the density profiles within dense cores. However, the observed α is sometimes steeper. For instance, molecular clouds like Polaris and Pipe have power-law indices α = 3.9 and α = 3.0, respectively (Lombardi et al. 2015). These clouds are diffuse and have much less star formation compared to active star-forming clouds like Aquila, which have α ≈ 2 (Könyves et al. 2015). Auddy et al. (2018) showed that the magnetic field can significantly affect the slope of the power-law tail. Clouds with a strong magnetic field (subcritical mass-to-flux ratio) and small amplitude initial perturbations develop a steep power-law tail (α ≈ 4), consistent with gravitationally driven ambipolar diffusion leading to shallower core density profiles than in a hydrodynamic collapse. In contrast, turbulent subcritical clouds retain the lognormal shape for a long time and eventually develop a power-law tail with α ≈ 2 in a region that has become supercritical due to turbulence-enhanced ambipolar diffusion.

Each PDF has at least three measurable parameters: the width of the lognormal part, the slope of the power-law tail, and the transitional column density Σ_TP that separates the lognormal from the power-law portion. Many theoretical studies have associated the standard deviation σ of the lognormal distribution with the sonic Mach number of driven turbulence (e.g., Federrath et al. 2008; Molina et al. 2012). The power-law tail develops when self-gravity is introduced into a driven turbulence simulation (Collins et al. 2012; Federrath & Klessen 2013). Decaying hydrodynamic turbulence simulations with self-gravity also show a rapid development of a power-law tail in the PDF (Ballesteros-Paredes et al. 2011; Kritsuk et al. 2011; Ward et al. 2014). Previous studies have not developed a theory for the location of Σ_TP, although Burkhart et al. (2017) and Imara & Burkhart (2016) have proposed that it is associated with the H i-to-H₂ transition in the interstellar medium. This has some appeal since H i clouds are known to be non-self-gravitating, whereas molecular clouds exist at higher pressures and are considered to be self-gravitating (Blitz 1991). However, observations of many molecular clouds show that the transition occurs within the molecular gas, and that the value of Σ_TP is unique to each cloud, e.g., NGC 3603, Carina, Maddalena, and Auriga all have different deviation points as listed in Table 1 of Schneider et al. (2015b). This implies that Σ_TP is an imprint of initial conditions inherent to a particular star-forming cloud and is set by physical processes. Furthermore, the longstanding well-known low efficiency of star formation (e.g., Goldsmith et al. 2008) within molecular clouds means that most of the molecular gas mass is not within gravitationally contracting dense cores that account for the power-law portion of the PDF.

In this Letter, we focus on the physical origin of the transition point Σ_TP, using a different approach than adopted in most previous studies of the PDF. Auddy et al. (2018) showed that in a decaying turbulence scenario with supercritical mass-to-flux ratio, the lognormal body is quickly lost and a power law is developed for essentially all densities past the peak. Simulations with constant turbulent driving with Fourier space perturbations are able to maintain a distinct lognormal body and a power-law tail in the PDF, although that may be an artifact of turning on self-gravity only after a steady turbulent driving has been established. Since real molecular clouds do not have a "switch-on" gravity, we investigate here the scenario that a strong magnetic field (i.e., subcritical mass-to-flux ratio) supports large amplitude oscillations (even while there is an overall decay of turbulence) that result in the maintenance of a lognormal-like body of the PDF. Gravity is always at work, but can only win out in dense regions that have undergone a rapid turbulence-accelerated ambipolar diffusion. In this case, Auddy et al. (2018) showed that a power-law tail with α ≈ 2 is added to the lognormal-like body of the PDF. Here, we develop a self-consistent theory of the origin of Σ_TP in this scenario, and test it against a suite of simulations with different initial conditions. We find a direct link between Σ_TP and the relative importance of turbulence and magnetic fields in the initial cloud. In Section 2 we present numerical simulations that study the properties of the column density PDFs. In Section 3 we derive an analytic expression for Σ_TP based on a model of turbulent compression of a magnetized cloud and compare with the numerical results. In Sections 4 and 5 we discuss and summarize our results, respectively.

We study the time evolution of the column density PDFs for five different models using three-dimensional magnetohydrodynamic (MHD) simulations including self-gravity and ambipolar diffusion. The numerical setup is similar to the ones previously used in Kudoh et al. (2007), Kudoh & Basu (2011), and Auddy et al. (2018), and we run with a number of grid points in each direction (N_x, N_y, N_z) = (512, 512, 20).

We consider models with a subcritical initial mass-to-flux ratio so that the magnetic field strength is dynamically important. The initial turbulent flow field causes the PDF to have a predominantly lognormal shape. However, as the cloud evolves, it forms compressed regions due to the large-scale flow and develops pockets of high column density. It then rebounds and shows oscillations. With each successive compression, more regions with high column density develop and cause a gradual widening of the width of the lognormal, even though the turbulence is decaying (see also Tassis et al. 2010; Ward et al. 2014; Auddy et al. 2018).

We follow the time evolution of the column density PDF for each model with different initial conditions. The PDF evolves over time from primarily having a lognormal shape at early times to developing a power-law tail (final time $\sim 1\,\mathrm{Myr}$) when gravity dominates. The power-law develops after several oscillations, as the local pockets of higher column density become supercritical and go into a runaway collapse.

2.1. Numerical Parameters

The initial state has a uniform density in x, y and is stratified in the z direction with a scale length ${H}_{0}={c}_{{\rm{s}}0}/\sqrt{2\pi G{\rho }_{0}}$, where c_s0 and ρ₀ are the isothermal sound speed and density at the midplane z = 0. For more details of the initial setup, see Kudoh et al. (2007), Kudoh & Basu (2011), and Auddy et al. (2018). We choose H₀, c_s0, and ρ₀ as units of length, velocity, and density, respectively. This gives the unit of time ${t}_{0}\equiv {H}_{0}/{c}_{{\rm{s}}0}$. The ratio of the initial gas to magnetic pressure at z = 0 is

$$\begin{eqnarray}&&{\beta }_{0}=\displaystyle \frac{8\pi {c}_{{\rm{s}}0}^{2}{\rho }_{0}}{{B}_{0}^{2}},\end{eqnarray} \tag{ 1 }$$

where B₀ is the initially uniform vertical magnetic field. We input Gaussian random velocity fluctuations of amplitude v_a for each of the x and y components of velocity, in which the Fourier spectrum is ${v}_{k}^{2}\propto {k}^{-4}$. Appropriate choices of ρ₀ and c_s0 lead to dimensional values of standard quantities. For example, if ${n}_{0}\equiv {\rho }_{0}/{m}_{n}={10}^{4}$ cm⁻³, where m_n = 2.33 × 1.67 × 10⁻²⁴ g and c_s0 = 0.2 km s⁻¹, we get H₀ ≃ 0.05 pc and t₀ ≃ 2.5 × 10⁵ yr. If β₀ = 0.16, Equation (1) yields B₀ ≃ 50 μG. The initial column density is Σ₀ = ρ₀H₀ ≃ 6 × 10⁻³ g cm⁻²; therefore, the number column density is ${N}_{0}\equiv {{\rm{\Sigma }}}_{0}/{m}_{n}\simeq 1.5\times {10}^{21}$ cm⁻².

2.2. Fitting Functions

In order to characterize the shape of the PDFs including the transition from the lognormal to power-law tail we consider two fitting functions: a purely lognormal function and a piecewise function that is a combination of a lognormal and a power law. If f(η) is the PDF, the lognormal model is

$$\begin{eqnarray}&&f{(\eta )}_{\mathrm{LN}}=\mathrm{log}\left[A\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\exp \left(-\displaystyle \frac{{(2.3\eta -\mu )}^{2}}{2{\sigma }^{2}}\right)\right],\end{eqnarray} \tag{ 2 }$$

where $\eta =\mathrm{log}({\rm{\Sigma }}/{{\rm{\Sigma }}}_{0})$, $A=\mathrm{ln}(10)\times {N}_{\mathrm{total}}\times {\rm{\Delta }}\mathrm{log}({\rm{\Sigma }}/{{\rm{\Sigma }}}_{0})$ is⁴ the normalization constant, μ is the mean, and σ is the standard deviation. The data are binned with a uniform spacing of ${\rm{\Delta }}\mathrm{log}({\rm{\Sigma }}/{{\rm{\Sigma }}}_{0})\simeq 0.02$. For the piecewise function (see also Myers 2015; Pokhrel et al. 2016) we consider a combination of a lognormal and a power law:

$$\begin{eqnarray}\begin{array}{rcl}f{(\eta )}_{\mathrm{LNPL}} & = & f{(\eta )}_{\mathrm{LN}},\quad \mathrm{if}\,\eta \leqslant {\eta }_{\mathrm{TP}},\\ & = & f{({\eta }_{\mathrm{TP}})}_{\mathrm{LN}}+\alpha \eta ,\quad \mathrm{if}\,\eta \gt {\eta }_{\mathrm{TP}},\end{array}\end{eqnarray} \tag{ 3 }$$

where α is the index of the power law and ${\eta }_{\mathrm{TP}}\equiv \mathrm{log}({{\rm{\Sigma }}}_{\mathrm{TP}}/{{\rm{\Sigma }}}_{0)})$ is the logarithmic value of the transition column density. Thus, we have two fitting functions and four free parameters: μ, σ, η_TP, and α.

We use this piecewise four-parameter function rather than a mathematically simpler three-parameter continuous function, for example, the modified lognormal power-law distribution (Basu et al. 2015), since it clearly identifies a transition point η_TP. We are then also using the same means to identify the transition point as used in observational analyses like Pokhrel et al. (2016).

2.3. Fitting the Simulation Data

The value of the fit parameters is essential to characterize the shape of the column density PDFs. However, the column density PDFs are evolving in time with their shape changing from purely lognormal to a hybrid function. Thus, it is essential to have a robust fitter that can capture the transition and identify η_TP. We fit the column density PDFs at different times using f(η)_LN (Equation (2)) and f(η)_LNPL (Equation (3)) and compute the resulting χ² values. We accept f(η)_LNPL only when its χ² value is less than 20% of that of the simpler lognormal function f(η)_LN and the power-law index α < 5. For models that are fit with f(η)_LNPL, we further use a Markov Chain Monte Carlo (MCMC) method (van Dyk 2003) as a second fitter. This gives us more robust best-fitting values from the parameter space along with reliable uncertainties. We have used the Python package Pymc for this purpose (Patil et al. 2010).

The power-law tail appears during the final stages of the simulations, primarily due to ambipolar-diffusion-driven gravitational contraction. All the free parameters evolve moderately, including the transitional column density η_TP, which grows by ≈10% from its initial appearance until the final time step. For simplicity we only consider the column density PDFs at the end of the simulation when the maximum density has reached 100ρ₀. Runaway collapse has ensued in the high-density regions at this time. While we cannot follow the PDFs into the protostellar phase in these simulations, we anticipate that the large-scale maps of the PDFs will remain largely the same.

Figure 1 shows the column density PDFs of five models with different initial conditions (i.e., β₀ and v_t0) along with the best-fit lognormal and power-law functions. The results are summarized in Table 1. The run time for each simulation is indicated on the top right of each plot. The best-fitting parameters α and η_TP are also shown. The plots are arranged from top left to bottom right according to increasing values of ${{ \mathcal M }}_{0}^{2}{\beta }_{0}$. The black dotted line marks the column density value at which the lognormal PDF ends and the power-law tail begins. The green shaded region shows the standard deviation of the η_TP value obtained from the MCMC fit.

Zoom In Zoom Out Reset image size

Table 1.  Model and Fit Parameters

Model	va/cs	β0	${{ \mathcal M }}_{0}^{2}{\beta }_{0}$	ηTP	$| \alpha |$	μ	σ
T1	2.0	0.16	1.28	0.62 ± 0.02	2.97 ± 0.10	0.59 ± 0.01	0.31 ± 0.01
T2	2.0	0.25	2.00	0.73 ± 0.03	2.68 ± 0.15	0.64 ± 0.02	0.35 ± 0.01
T3	3.0	0.16	2.88	0.79 ± 0.01	2.72 ± 0.09	0.55 ± 0.01	0.41 ± 0.01
T4	2.5	0.36	4.50	1.01 ± 0.03	1.78 ± 0.32	0.63 ± 0.01	0.46 ± 0.01
T5	3.0	0.36	6.48	1.20 ± 0.02	1.64 ± 0.31	0.53 ± 0.01	0.57 ± 0.01

Note. Fit parameters for the piecewise lognormal and power-law function are α, σ, μ, and η_TP. β₀ is the initial ratio of thermal to magnetic pressure at z = 0, v_a is the amplitude of the initial velocity fluctuation, and ${{ \mathcal M }}_{0}$ is the sonic Mach number.

Download table as:  ASCIITypeset image

The transition point η_TP shifts toward higher column density with an increasing value of ${{ \mathcal M }}_{0}^{2}{\beta }_{0}$. For example, it is minimum for model T1 (${{ \mathcal M }}_{0}^{2}{\beta }_{0}=1.28$) and maximum for model T5 (${{ \mathcal M }}_{0}^{2}{\beta }_{0}=6.48$) with η_TP = 0.62 ± 0.02 and η_TP = 1.20 ± 0.02, respectively. The value depends on the strength of the initial magnetic field and the amplitude of the velocity perturbation.

We can understand the physical origin of the transition point with an analytic model in which the magnetic field is dynamically important. The cloud flattens along the mean magnetic field direction (z), and the subsequent evolution is primarily perpendicular to the magnetic field. The mass-to-flux ratio is subcritical until ambipolar diffusion creates supercritical pockets that are prone to collapse. Turbulence causes the creation of locally compressed regions that have a pressure balance between magnetic and ram pressure. This results in the formation of magnetic ribbons (Auddy et al. 2016; see also Kudoh & Basu 2014). The cloud is stratified along the z direction with compression along the x–y plane. We simplify the analysis and assume that the thermal pressure is negligible compared to the magnetic pressure and the ram pressure of the flow. The pressure due to the magnetic field B upon compression balances the initial pressure due to the background magnetic field B₀ and the external ram pressure in the x–y direction:

$$\begin{eqnarray}&&H\displaystyle \frac{{B}^{2}}{8\pi }={H}_{0}\left({\rho }_{0}{v}_{{\rm{t}}0}^{2}+\displaystyle \frac{{B}_{0}^{2}}{8\pi }\right),\end{eqnarray} \tag{ 4 }$$

where ${v}_{{\rm{t}}0}=\sqrt{2}{v}_{a}$ is the nonlinear flow speed. This results in a quasi-equilibrium state as compression ceases and oscillations begin. The gas has already settled into a hydrostatic equilibrium along the z direction and the cloud has a half-thickness

$$\begin{eqnarray}&&H=\displaystyle \frac{{c}_{{\rm{s}}}}{\sqrt{2\pi G\rho }}\end{eqnarray} \tag{ 5 }$$

(Spitzer 1942). Integrating the density along the scale height H in the z direction gives the column density

$$\begin{eqnarray}&&{\rm{\Sigma }}=2\rho H.\end{eqnarray} \tag{ 6 }$$

The initial density compression is very high (see Figure 1 in Kudoh & Basu 2008), but the subcritical mass-to-flux ratio results in a strong rebound. We consider that the cloud is nearly flux frozen during its initial compression, i.e., B/Σ = constant, as the ambipolar diffusion time is much longer than the compression time. Using (5) and (6) in Equation (4) along with the flux-frozen condition we find

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}=\left[{v}_{{\rm{t}}0}^{2}\left(\displaystyle \frac{8\pi {\rho }_{0}}{{B}_{0}^{2}}\right)+1\right].\end{eqnarray} \tag{ 7 }$$

The force balance of Equation (7) gives a critical column density that we denote as Σ_TP. However, subsequent oscillations are not as strong due to decay of the initial turbulence amplitude v_t0 and loss of magnetic flux due to ambipolar diffusion. Allowing for such variations we rewrite Equation (7) in terms of the sonic Mach number ${{ \mathcal M }}_{0}={v}_{{\rm{t}}0}/{c}_{{\rm{s}}}$, and plasma β₀ as

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{TP}}}{{{\rm{\Sigma }}}_{0}}=a({{ \mathcal M }}_{0}^{2}{\beta }_{0}+1),\end{eqnarray} \tag{ 8 }$$

where a is a correction factor of order unity that contains uncertainties about the flux loss and turbulent decay.

Here ${{ \mathcal M }}_{0}^{2}{\beta }_{0}\equiv 2{{ \mathcal M }}_{{\rm{A}}0}^{2}$, where ${{ \mathcal M }}_{{\rm{A}}0}={v}_{{\rm{t}}0}/{v}_{{\rm{A}}0}$ is the initial Alfvénic Mach number and ${v}_{{\rm{A}}0}={B}_{0}/\sqrt{(4\pi {\rho }_{0})}$ is the initial Alfvén speed in the midplane.

3.1. Physical Interpretation of Σ_TP

Figure 2 shows the variation of the normalized transitional column density (Σ_TP/Σ₀) for simulated models with different initial values of ${{ \mathcal M }}_{0}^{2}{\beta }_{0}$. We fit the analytic expression (Equation (8)) to the simulation data and get a best-fit value of a = 1.9. We find a good agreement between the simulation data and our analytic model.

Zoom In Zoom Out Reset image size

The transition from the lognormal to the power-law tail signifies both structural and morphological changes. It marks a transition from the ambient subcritical turbulent background (Σ < Σ_TP) to a compressed denser region where Σ > Σ_TP. Due to ambipolar diffusion (ion–neutral drift) the force balance between the ram pressure and the magnetic field gradually relaxes. There is a gradual loss of magnetic flux as the neutrals diffuse past the ions with each successive oscillation. The density is enhanced after each compression resulting in an increase of mass-to-flux ratio. The transitional column density Σ_TP defines this cutoff beyond which the mass-to-flux ratio becomes critical. For Σ > Σ_TP gravity becomes increasingly important and the power-law tail emerges.

Furthermore, with increased strength of the Alfvénic Mach number ${{ \mathcal M }}_{{\rm{A}}0}$ the initial compression is much stronger and it results in a higher density. This causes Σ_TP to shift toward higher values with increasing strength of $2{{ \mathcal M }}_{{\rm{A}}0}^{2}$.

The shape of the column density PDF is an imprint of the initial conditions in a star-forming molecular cloud. It holds the key to finding the link between the structural properties and the ambient conditions that trigger star formation in molecular clouds. However, it is difficult to detect low column density material (the lognormal part) using dust emission and extinction measurement because of observational biases (Lombardi et al. 2015). The lognormal peak is sometimes considered an artifact arising due to data incompleteness (Alves et al. 2017) or undetectable due to insufficient sampling or limited field of view (Körtgen et al. 2019). The more robust observational quantities are the characteristic break Σ_TP in the PDF and the power-law slope (α), as these are less affected by such constraints (Lombardi et al. 2015).

Our model shows that Σ_TP is a measure of the initial Alfvénic Mach number ${{ \mathcal M }}_{{\rm{A}}0}$. Observations of Σ_TP can be used to infer the value of ${{ \mathcal M }}_{{\rm{A}}0}\propto {{ \mathcal M }}_{0}{\beta }_{0}^{1/2}$. This could in principle also lead to an estimate of the initial normalized mass-to-flux ratio μ₀, where ${\mu }_{0}^{2}\simeq {\beta }_{0}$ (see Kudoh et al. 2007), if the initial Mach number ${{ \mathcal M }}_{0}$ can be estimated. This is potentially important since a direct measurement of the magnetic field strength using the Zeeman effect is difficult (Crutcher 2012). Indirect probes of the magnetic field such as dust polarization (Hoang & Lazarian 2008), spectroscopic methods (Auddy et al. 2019), and Faraday rotation (Wolleben & Reich 2004) also have their limitations. An estimate of ${{ \mathcal M }}_{0}$ can in principle be made from the observed width σ of the lognormal PDF by developing an analytic/empirical relation that captures the increase of σ with the increasing strength of ${{ \mathcal M }}_{0}^{2}{\beta }_{0}$ (see Table 1). Such an analysis can be pursued in future work.

Our key findings are:

• 1.  
  The transitional column density Σ_TP represents a transition from a turbulent magnetically dominated background (Σ < Σ_TP) having lognormal shape to a dense region (Σ > Σ_TP) with a power-law tail where gravity is dominant.
• 2.  
  Σ_TP marks the boundary between regions with subcritical (magnetically dominated) and supercritical (weak magnetic field) mass-to-flux ratio in a star-forming molecular cloud.
• 3.  
  Σ_TP depends on the initial velocity perturbation (sonic Mach number ${{ \mathcal M }}_{0}$) and the magnetic field strength (plasma β₀). Alternatively, it is a measure of the initial Alfvénic Mach number and increases with the increasing strength of $2{{ \mathcal M }}_{{\rm{A}}0}^{2}$.

We thank the anonymous referee for constructive comments. S.A. acknowledges support from an ASIAA Postdoctoral Fellowship. Computations were carried out using facilities of SHARCNET. S.B. is supported by a Discovery Grant from NSERC.