The HI-to-H2 Transition in a Turbulent Medium

TL;DR

This study investigates how turbulence-induced density fluctuations affect the atomic-to-molecular hydrogen transition in photodissociation regions, using analytical models and numerical simulations to interpret HI observations.

Contribution

It introduces an analytic model for the HI PDF in turbulent media and demonstrates its application to observational data, linking turbulence properties to HI structure.

Findings

Mean HI column density remains consistent with uniform-density models.

HI PDF width depends on radiation, Mach number, and turbulence scale.

Application to Perseus cloud suggests small-scale turbulence influences HI structure.

Abstract

We study the effect of density fluctuations induced by turbulence on the HI/H$_2$ structure in photodissociation regions (PDRs) both analytically and numerically. We perform magnetohydrodynamic numerical simulations for both subsonic and supersonic turbulent gas, and chemical HI/H$_2$ balance calculations. We derive atomic-to-molecular density profiles and the HI column density probability density function (PDF) assuming chemical equilibrium. We find that while the HI/H$_2$ density profiles are strongly perturbed in turbulent gas, the mean HI column density is well approximated by the uniform-density analytic formula of Sternberg et al. (2014). The PDF width depends on (a) the radiation intensity to mean density ratio, (b) the sonic Mach number and (c) the turbulence decorrelation scale, or driving scale. We derive an analytic model for the HI PDF and demonstrate how our model, combined with 21 cm observations, can be used to constrain the Mach number and driving scale of turbulent gas. As an example, we apply our model to observations of HI in the Perseus molecular cloud. We show that a narrow observed HI PDF may imply small scale decorrelation, pointing to the potential importance of subcloud-scale turbulence driving.

Figures

Equations

$$Rn{n}_{\mathrm{HI}}=\frac{1}{2}{D}_0{f}_{\mathrm{shield}}(N_{{\mathrm{H}}_2}){\mathrm{e}}^{-{\sigma }_{g}N}{n}_{{\mathrm{H}}_2},$$

$$\displaystyle N_{{\rm HI}}\$$

$$=5.3\times 10^{20}\frac{1}{{\varphi }_g{Z}^{\prime }}\ln \left[\frac{\alpha G}{2}+1\right]\mathrm{c}{\mathrm{m}}^{-2}$$

$${\sigma }_g=1.9\times 10^{-21}{\varphi }_g{Z}^{\prime }\mathrm{c}{\mathrm{m}}^2$$

$$\alpha G=\frac{D_{0}G}{Rn}=2.0I_{\mathrm{UV}}\left(\frac{30\mathrm{c}{\mathrm{m}}^{-3}}{n}\right),$$

$$t_{\mathrm{chem}}=\frac{1}{2Rn+D}.$$

$$t_{\mathrm{chem}}=\frac{1}{4Rn}\approx 2.7\frac{1}{{\varphi }_R{Z}^{\prime }{T}_2^{1/2}{n}_2}\mathrm{Myr},$$

$$R=3\times 10^{-17}{T}_2^{1/2}{\varphi }_R{Z}^{\prime }\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1},$$

$$n\ge \frac{1}{4R{t}_{\mathrm{cloud}}}\approx 27\frac{1}{{\varphi }_R{Z}^{\prime }{T}_2^{1/2}}\left(\frac{10\mathrm{Myr}}{t_{\mathrm{cloud}}}\right)\mathrm{c}{\mathrm{m}}^{-3},$$

$$t_{\mathrm{turb}}\simeq \frac{L_{\mathrm{HI}}}{\delta {\text{varv}}_{1d}(L_{\mathrm{HI}})}$$

$$L_{\mathrm{HI}}=\frac{1}{{\sigma }_{g}n}=1.7\frac{1}{{\varphi }_g{Z}^{\prime }{n}_2}\mathrm{pc},$$

$$\delta \text{varv}(\ell )=\delta \text{varv}(L_{\mathrm{drive}}){\left(\frac{\ell }{L_{\mathrm{drive}}}\right)}^{1/2}(L_s\le \ell \le L_{\mathrm{drive}}),$$

$$t_{\mathrm{turb}}$$

$$\approx 4.0\frac{1}{\text{pazocal}{M}_s{T}_2^{1/2}{n}_2{\varphi }_g{Z}^{\prime }}{\left(\frac{L_{\mathrm{HI}}}{L_{\mathrm{drive}}}\right)}^{-1/2}\mathrm{Myr},$$

$$\frac{t_{\mathrm{chem}}}{t_{\mathrm{turb}}}\approx 0.66\text{pazocal}{M}_s\frac{{\varphi }_g}{{\varphi }_R}{\left(\frac{L_{\mathrm{HI}}}{L_{\mathrm{drive}}}\right)}^{1/2}.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \text{varv})$$

$$\frac{\partial \rho \text{varv}}{\partial t}+\nabla \cdot \left[\rho \text{varv}\text{varv}+\left(p+\frac{B^2}{8\pi }\right)\mathbf{I}-\frac{1}{4\pi }\mathbf{B}\mathbf{B}\right]$$

$$\frac{\partial \mathbf{B}}{\partial t}-\nabla \times (\text{varv}\times \mathbf{B})$$

$${\sigma }_x$$

$${\sigma }_{\ln x}^2$$

$$x_{\ell }\equiv \frac{{\int }_0^{\ell }x\mathrm{d}{\ell }^{\prime }}{\ell },$$

$$\text{pazocal}N(\ell )=1+\frac{\ell }{L_{\mathrm{dec}}},$$

$${\sigma }_{x_{\ell }}\simeq \frac{{\sigma }_x}{\sqrt{\text{pazocal}N(\ell )}}.$$

$$\frac{{\sigma }_{x_l}}{{\sigma }_x}={\left(1+\frac{\ell }{L_{\mathrm{dec}}}\right)}^{-1/2}={\left(1+\frac{L_{\mathrm{box}}}{L_{\mathrm{dec}}}\frac{\ell }{L_{\mathrm{box}}}\right)}^{-1/2},$$

$$\frac{L_{\mathrm{dec}}}{L_{\mathrm{box}}}=7.8\times 10^{-2},$$

$$\frac{L_{\mathrm{dec}}}{L_{\mathrm{drive}}}=0.20.$$

$$\overline{\alpha G}\equiv \frac{D_{0}G}{R⟨n⟩}=2.0I_{\mathrm{UV}}\left(\frac{30\mathrm{c}{\mathrm{m}}^{-3}}{⟨n⟩}\right),$$

$$L_{\mathrm{HI}}\equiv \frac{1}{{\sigma }_g⟨n⟩}=5.7\left(\frac{30\mathrm{c}{\mathrm{m}}^{-3}}{⟨n⟩}\right)\frac{1}{{\varphi }_g{Z}^{\prime }}\mathrm{pc}.$$

$$\frac{L_{\mathrm{dec}}}{L_{\mathrm{HI}}}={\sigma }_g⟨n⟩L_{\mathrm{dec}}\equiv {\tau }_{\mathrm{dec}},$$

$$\text{pazocal}N(L_{\mathrm{HI}})=1+\frac{L_{\mathrm{HI}}}{L_{\mathrm{dec}}}=1+{\tau }_{\mathrm{dec}}^{-1},$$

$$L_{\mathrm{drive}}=5.0{\tau }_{\mathrm{dec}}{L}_{\mathrm{HI}}.$$

$$N_{\mathrm{HI}}=\frac{1}{{\sigma }_g}\ln \left[\frac{\overline{\alpha G}}{2}\frac{1}{x_L}+1\right].$$

$$\frac{\mathrm{d}f}{\mathrm{d}\ln N_{\mathrm{HI}}}=\left(1+\frac{2x_L}{\overline{\alpha G}}\right)\ln \left(\frac{\overline{\alpha G}}{2x_L}+1\right)\frac{\mathrm{d}f}{\mathrm{d}\ln x_L},$$

Full text

††slugcomment: Submitted to ApJ, 20 March 2017

1. Introduction

Giant molecular clouds serve as the nurseries for new stars in our Galaxy and in external galaxies (McKee & Ostriker, 2007). On global scales, observations of CO and dust show that the star-formation rate (SFR) surface density ($\Sigma_{\rm SFR}$) correlates with the H2 mass surface density ($\Sigma_{\rm H_{2}}$), following an almost linear trend (Bigiel et al., 2008; Genzel et al., 2010; Schruba et al., 2011; Tacconi et al., 2013; Azeez et al., 2016). The presence of H2 molecules is a basic ingredient for the formation of other heavy molecules such as CO, OH and H2O that serve as efficient coolants of cold gas (e.g., Herbst & Klemperer, 1973; Sternberg & Dalgarno, 1995; Tielens, 2013; van Dishoeck et al., 2013; Bialy & Sternberg, 2015). The study of far-ultraviolet (UV) shielding and the subsequent HI-to-H2 conversion is of fundamental importance for star-formation and molecule formation in the interstellar medium (ISM).

The HI-to-H2 transition in the interstellar medium of galaxies has been investigated by numerous authors over the last several decades, through analytic and numerical modeling (e.g., Federman et al., 1979; van Dishoeck & Black, 1986; Sternberg, 1988; Kaufman et al., 1999; Goldsmith et al., 2007; Gnedin & Draine, 2014; Liszt, 2015), as well as via hydrodynamics simulations (e.g., Robertson & Kravtsov, 2008; Gnedin et al., 2009; Glover et al., 2010; Dave et al., 2013; Thompson et al., 2014; Lagos et al., 2015; Hu et al., 2016) and observations (e.g., Savage et al., 1977; Reach et al., 1994; Rachford et al., 2002; Gillmon & Shull, 2006; Lee et al., 2012; Balashev et al., 2014; Noterdaeme et al., 2015; Nakanishi & Sofue, 2016). Analytic treatments of the HI-to-H2 transition have been presented by Krumholz et al. (2008), McKee & Krumholz (2010) and Sternberg et al. (2014) using a Strömgren type analysis for the total steady state column density of HI that is maintained by an incident photodissociating flux. In particular, Sternberg et al. (2014, hereafter, S14) derived a scaling relationship for the total HI column density in optically thick uniformly irradiated slabs as a function of the far-UV flux, the gas density, the dust absorption cross section and the H2 formation rate coefficient. Bialy & Sternberg (2016, hereafter BS16) presented an analytic procedure for generating atomic (${\rm HI}$) to molecular (${\rm H_{2}}$) density profiles for optically thick hydrogen gas clouds in Galactic star-forming regions. These studies thus far have been instrumental in interpreting emission line observations of HI /H2 interfaces (Lee et al., 2012; Bialy et al., 2015; Burkhart et al., 2015; Bihr et al., 2015; Bialy et al., 2017; Maier et al., 2017), for estimating star-formation thresholds in external galaxies (Leroy et al., 2008; Lada et al., 2012; Clark & Glover, 2014; Bialy & Sternberg, 2016; Burkhart & Loeb, 2016), and for sub-grid components in hydrodynamics simulations (Pelupessy et al., 2006; Thompson et al., 2014; Tomassetti et al., 2015).

Despite the progress towards an analytic theory for the physics of the HI-to-H2 transition, no current theory includes realistic turbulent density fluctuations. The turbulent nature of molecular clouds is evident from a variety of observations including non-thermal broadening (Stutzki & Guesten, 1990; Dickey et al., 2001; Heiles & Troland, 2003; Heyer & Brunt, 2004), velocity/density power spectrum (Stanimirović & Lazarian, 2001; Swift & Welch, 2008; Chepurnov & Lazarian, 2009; Pingel et al., 2013; Chepurnov et al., 2015), and fractal and hierarchical structures (Elmegreen & Elmegreen, 1983; Vazquez-Semadeni, 1994; Burkhart et al., 2013a). Simulations and observations have shown that supersonic turbulence creates filaments and regions of high density contrast (Kowal et al., 2007; Burkhart et al., 2009; Federrath et al., 2010). This behavior suggests that the assumption of uniform density in current analytic models for the HI-to-H2 transition should be revisited.

The effects of turbulence on the chemical structure of interstellar clouds has been studied from various perspectives. Xie et al. (1995), Willacy et al. (2002) and Bell et al. (2011) studied the effects of turbulent mixing of chemical species through a diffusion approximation. They found that atomic abundances (e.g., H, C, and O) may be significantly increased in the interiors of molecular clouds if the diffusion coefficient is large. Levrier et al. (2012) studied the chemical structure of turbulent photodissociation regions (PDRs) using a post-processing approach, and found that the abundances of various molecules (e.g., H2, CO, CH, and CN) exhibit strong deviations from a homogeneous PDR model (cf., Offner et al., 2013). Glover & Mac Low (2007), Glover et al. (2010), Micic et al. (2011) and Valdivia et al. (2016) performed MHD simulations and followed the time-dependent H2 formation self-consistently. They focused on the molecular content and showed that strong density compressions created by supersonic turbulent flows produce H2 rapidly on timescales of few Myrs.

In this paper we study the effects of turbulent density perturbations on the ${\rm HI}/{\rm H_{2}}$ structure in PDRs, focusing on the atomic gas produced by photodissociation at the cloud boundaries. As we show, this gas is particularly useful for constraining the nature of turbulence, via 21 cm observations. We consider a twofold approach, (i) via numerical MHD simulations supplemented by ${\rm H_{2}}/{\rm HI}$ chemical balance calculation, and (ii) analytic modeling, and introducing a novel method for constraining the Mach number and the turbulence driving scale.

The paper is organized as follows: In §2 we provide a basic overview of HI-to-H2 theory for uniform-density gas, as presented by S14 and BS16. In §3 we discuss the effect of density fluctuations and the validity of our chemical steady-state assumption when considering turbulence. In §4 we present the results of our MHD simulations. In §5 we present computations of ${\rm HI}/{\rm H_{2}}$ profiles and HI column density PDFs for turbulent media. In §6 we develop an analytic model for the HI column density distribution. In §7 we demonstrate how our analytic model may be used to constrain turbulent parameters from 21 cm observations. We discuss and summarize our results in §8.

2. Uniform density gas

In this section we review briefly the theory of HI-to-H2 transition in steady-state, uniform density gas. For a through discussion we refer the reader to S14 and BS16.

At any cloud depth, and for unidirectional radiation normal to the cloud surface, H2 formation-destruction equilibrium is given by

[TABLE]

where $R$ (cm3 s*-1*) is the H2 formation rate coefficient, $n=n_{{\rm HI}}+2n_{{\rm H_{2}}}$ is the total (atomic plus molecular) hydrogen volume density and $D_{0}$ (s*-1*) is the free space H2 photodissociation rate. In this expression, $f_{\rm shield}$ is the H2 self-shielding function that depends on the H2 column density $N_{{\rm H_{2}}}$ and also on the absorption-line Doppler broadening parameter $b_{D}$ (km s*-1*) (Draine & Bertoldi 1996, S14). The factor $\mathrm{e}^{-\sigma_{g}N}$ is the dust absorption attenuation term, where $\sigma_{g}$ (cm2) is the dust-grain absorption cross section per hydrogen nuclei integrated over the Lyman-Werner dissociation band (11.2-13.6 eV; hereafter LW-band), and $N=N_{{\rm HI}}+2N_{{\rm H_{2}}}$ is the total, atomic plus molecular column density. The factor of 1/2 accounts for absorption of half the radiation by the optically thick slab.

Assuming that all of the photodissociating radiation is absorbed in the cloud, the total HI column density, converges to a finite value, $N_{{\rm HI}}$. As shown by S14, assuming slab geometry, and for constant density

[TABLE]

where in the second equality

[TABLE]

where $Z^{\prime}$ is the dust-to-gas ratio relative to Galactic, and $\phi_{g}$ is a factor of order unity that characterizes the dust absorption properties. In Eq. (2), $\alpha\equiv D_{0}/(Rn)$ is the (dimensionless) ratio of free-space H2 photodissociation and H2 formation rates. The dimensionless factor $G\equiv\sigma_{g}\int f_{\rm shield}(N_{{\rm H_{2}}})\mathrm{e}^{-2\sigma_{g}N_{{\rm H_{2}}}}dN_{{\rm H_{2}}}$ is “the effective shielding factor” and may be expressed as $G=3.0\times 10^{-5}\phi_{g}Z^{\prime}(9.9/[1+8.9\phi_{g}Z^{\prime}])^{0.37}$ (BS16). The combination

[TABLE]

has the physical meaning of an effective dissociation parameter taking into account H2-shielding and the competition with dust absorption. The numerical value in Eq. (4) is for $R=3\times 10^{-17}$ cm3 s*-1*, $\phi_{g}Z^{\prime}=1$, and $D_{0}=5.8\times 10^{-11}\ I_{\rm UV}$ s*-1* (S14) where $I_{\rm UV}$ is the radiation strength relative to the Draine (1978) field. The product $\alpha G$ may be small or large for realistic astronomical environments. For example, for the star-forming region W43 and for the Perseus molecular cloud, Bialy et al. (2017) and Bialy et al. (2015) derived ${\alpha G}\sim 20$ and ${\alpha G}\sim 10$, respectively, whereas for a sample of dwarf irregular galaxies in the LITTLE THINGS survey (Hunter et al., 2012), Maier et al. (2017) deduced ${\alpha G}<1$. For cold neutral medium (CNM) which is in pressure equilibrium with the warm neutral medium (WNM), the $n/I_{\rm UV}$ ratio is restricted to the range $\approx 8-70$ cm*-3* (Wolfire et al., 2003), giving $({\alpha G})_{\rm CNM}\approx 1-8$.

The left panel of Fig. 1 shows the HI and ${\rm H_{2}}$ fractional abundance profiles, $x_{{\rm HI}}\equiv n_{{\rm HI}}/n$ and $2x_{{\rm H_{2}}}\equiv 2n_{{\rm H_{2}}}/n$, as functions of cloud depth, as parameterized by the gas column density $N$. The various curves are for $\phi_{g}Z^{\prime}=1$, ${\alpha G}=2$ and Doppler parameter ranging from $b_{D}=0.5$ to 8 km s*-1*. With increasing depth the radiation is absorbed by H2 photodissociation events and by dust absorption, and the gas makes the transition from atomic to molecular form. For larger $b_{D}$, the Doppler cores of the H2 lines are broader, and the onset of self shielding occurs at larger cloud depths. However, the HI-to-H2 transition point is insensitive to $b_{D}$ because it occurs deeper in the cloud where the LW-flux is absorbed in the H2 damping wings (BS16).

The right panel of Fig. 1 shows the total HI column density, $N_{{\rm HI}}$, as a function of ${\alpha G}$. For large ${\alpha G}$ (“the strong field limit”) dust absorption determines the HI column density, and $N_{{\rm HI}}$ is weakly dependent on ${\alpha G}$. For small $\alpha G$ (“the weak field limit”), H2 self-shielding dominates and $N_{{\rm HI}}\propto{\alpha G}$. As for the HI-to-H2 transition points, $N_{{\rm HI}}$ is insensitive to $b_{D}$ (and for the same reason). The insensitivity of $N_{{\rm HI}}$ to the Doppler parameter is important in our analysis of turbulent media and the line broadening induced by turbulent motions.

3. Density Fluctuations and Timescales

As discussed above, two basic assumptions for the HI-to-H2 density profiles and the HI column density (as shown in Fig. 1 and given by Eq. 2), are a constant gas density $n$ and chemical steady state. In this paper we relax the assumption of constant density by considering turbulent density fluctuations, while retaining the assumption of chemical steady state. Density fluctuations are naturally produced in the cold ISM by supersonic turbulence. Since the H2 formation-to-removal rate ratio is proportional to density, a local increase in density increases the local H2 fraction. Such a perturbation also affects deeper locations in the cloud through enhanced H2 self-shielding which depends on the H2 column density.

For a given density $n$, and a local (attenuated) dissociation rate $D$, the chemical time is

[TABLE]

In the outer atomic layers, $D\gg 2Rn$ and $t_{\rm chem}\simeq 1/D$ is the photodissociation time, which is typically very short. For example, for unshielded gas, $D=D_{0}/2=2.9\times 10^{-11}$ s*-1* and $t_{\rm chem}=1.1\times 10^{3}$ yr. Beyond the atomic-to-molecular transition, $D<2Rn$ and $t_{\rm chem}\simeq 1/(2Rn)$ is the H2 formation time, which can become long. The strongest effect of density perturbations will occur near the HI-to-H2 transition points where $Rn$ and $D$ are comparable. At the transition point, $D=2Rn$, and

[TABLE]

where the rate coefficient

[TABLE]

Here, $T_{2}\equiv T/10^{2}$ K, $n_{2}\equiv n/10^{2}$ cm*-3*, and $\phi_{R}$ is a factor of order unity. Typically, $\phi_{R}\approx 1$, however, in some environments the H2 formation rate may be enhanced. For example, Habart et al. (2003, 2004) found that moderately illuminated PDRs, such as Oph W, S140 and IC 63, may have $\phi_{R}\approx 5$, considerably reducing the chemical time.

Irrespective of any density fluctuations, the chemical time must be short compared to the cloud lifetime, $t_{\rm cloud}$. This requires

[TABLE]

where we have normalized $t_{\rm cloud}$ to characteristic lifetime of 10 Myrs.

In a turbulent medium, we also require $t_{\rm chem}\leq t_{\rm turb}$, where $t_{\rm turb}$ is the characteristic time over which turbulent density fluctuations are formed and destroyed. The turbulent time is

[TABLE]

where

[TABLE]

is the characteristic length-scale of the HI layer, and $\delta\varv_{1d}(L_{{\rm HI}})$ is the 1d velocity dispersion over $L_{{\rm HI}}$. Following the linewidth-size relation (Larson, 1981; McKee & Ostriker, 2007), the 3d velocity dispersion over a length-scale $\ell$ is

[TABLE]

where $L_{\rm drive}$ is the outer driving scale, and $L_{s}$ is the sonic length for which $\delta\varv(L_{s})\equiv c_{s}$, where $c_{s}$ is the sound speed. Defining the Mach number $\pazocal{M}_{s}\equiv\delta\varv(L_{\rm drive})/c_{s}$, and assuming an isotropic velocity field, we get

[TABLE]

where in the second equality we used $c_{s}=0.72T_{2}^{1/2}$ km s*-1*, assuming a mean particle mass of 1.6 the proton mass, at the transition point where $x_{{\rm HI}}=2x_{{\rm H_{2}}}$, and including Helium.

The ratio of the chemical and turbulent times is then

[TABLE]

Eqs. (12) and (13) are for $L_{s}\leq L_{{\rm HI}}\leq L_{\rm drive}$. For $L_{{\rm HI}}>L_{\rm drive}$, $\delta\varv\rightarrow\delta\varv(L_{\rm drive})$, and $L_{{\rm HI}}/L_{\rm drive}$ should be replaced with unity. For $L_{{\rm HI}}<L_{s}$ density perturbations are negligible as they are smoothed out by pressure waves. At $L_{{\rm HI}}=L_{s}$, $L_{{\rm HI}}/L_{\rm drive}=1/\pazocal{M}_{s}^{2}$. Thus, for moderate Mach numbers, $t_{\rm chem}/t_{\rm turb}$ remains close to unity, possibly ranging from $0.66(\phi_{g}/\phi_{R})$ to $0.66(\phi_{g}/\phi_{R})\pazocal{M}_{s}$. Throughout this paper we assume chemical equilibrium. We will consider the more complicated time-dependent problem elsewhere.

Turbulent motions may also affect the ${\rm HI}$/${\rm H_{2}}$ structure by shifting the frequencies of the H2 absorption lines, and reducing the efficiency of H2-self shielding. This affects the ${\rm HI}$/${\rm H_{2}}$ profiles at intermediate depths, $10^{14}\leq N_{{\rm H_{2}}}\leq 10^{18}$ cm*-2* where absorption is dominated by the Doppler cores (Gnedin & Draine, 2014). However, since most of the HI gas is accumulated at greater cloud depths where the radiation is absorbed in the H2-line damping wings, the velocity shifts do not affect $N_{{\rm HI}}$ (see the discussion in §2). The Doppler broadening is important for optically thin medium (to the LW-radiation), but not to optically thick clouds that have fully converted $\text{HI-to-H}_{2}$. Since our focus is on optically thick gas, we assume a constant $b_{D}=2$ km s*-1* throughout our calculations, and do not include any corrections to the Doppler parameter or the H2 self-shielding function (e.g., Gnedin & Draine, 2014).

4. MHD Simulations

In this Section, we use MHD simulations to obtain realistic density profiles for sub- and supersonic turbulent gas. Our simulations are isothermal and non-self gravitating (cf. Jappsen et al., 2005). This allows a natural extension of the S14 HI-to-H2 transition model, which is inherently isothermal, into the turbulent regime. Furthermore, an isothermal equation of state (EOS) allows a simple estimate of the density dispersion (Eq. 17, below). In reality, the HI-to-H2 transition takes place in a non-isothermal medium with heating and cooling processes acting, e.g., depth dependent photoelectric heating versus [CII] emission line cooling (Tielens & Hollenbach, 1985; Sternberg & Dalgarno, 1989). However, for moderate sonic Mach numbers ($\pazocal{M}_{s}\lesssim 5$), the density and column density PDFs, are similar in simulations of isothermal or non-isothermal EOS (e.g., Glover & Mac Low, 2007; Federrath & Banerjee, 2015).

We use a third-order-accurate hybrid essentially nonoscillatory scheme (Cho & Lazarian, 2002) to solve the ideal MHD equations,

[TABLE]

where $\rho$ is density, ${\bf B}$ is magnetic field, $p$ is the gas pressure, ${\bf I}$ is the identity matrix and $\bf{f}$ is the specific force. We assume zero-divergence condition $\nabla\cdot{\bf B}=0$, periodic boundary conditions, and an isothermal equation of state $p=c_{s}^{2}\rho$. For the source term $\bf{f}$, we assume a random large-scale solenoidal driving at a wave number $k\approx 2.5$ (i.e. 1/2.5 the box size). The simulations have 5123 resolution elements and have been employed in many previous works (Cho & Lazarian, 2003; Burkhart et al., 2009, 2010; Kowal et al., 2007, 2009, 2011).

Each simulation is defined by the sonic Mach number $\pazocal{M}_{s}\equiv|\boldsymbol{\varv}|/c_{s}$, and the Alfvénic Mach number $\pazocal{M}_{A}\equiv|\boldsymbol{\varv}|/\langle\varv_{A}\rangle$, where $\boldsymbol{\varv}$ is the velocity, $c_{s}$ and $\varv_{A}$ are the isothermal sound speed and the Alfvén speed, and $\langle\cdot\rangle$ denotes averages over the entire simulation box. We show results for $\pazocal{M}_{s}=0.5$, 2, and 4.5 simulations, i.e., subsonic, transonic, and supersonic gas. As we show below, the value of the sonic Mach number strongly affects the variance of the density field. The simulations are sub-alfvénic with $\pazocal{M}_{A}=0.7$ (i.e. strong magnetic field). We have also considered super-Alfvénic ($\pazocal{M}_{A}=2.0$) simulations and found that the results are weakly sensitive to the value of $\pazocal{M}_{A}$. Because the simulations are non self-gravitating they are scale-free and we may assign any desired physical scale for the box length and density (see Hill et al. 2008, Apendix). In this section we keep the results general and do not apply any physical scaling to the simulations. We scale the simulations to physical units in §5 below.

4.1. The 3D Density Distribution

Fig. 2 shows three random density cuts (upper panels) through the $\pazocal{M}_{s}=0.5,2$ and 4.5 simulations. The color axis corresponds to $\ln x$ where $x\equiv n/\langle n\rangle$. The density is nearly uniform for the subsonic simulation, but once the Mach number exceeds unity strong density fluctuations are generated. The lower panels show the PDFs of $\ln x$ for these simulations (shaded) . The $\ln x$ distributions are nearly Gaussian with a standard deviation that increases with Mach number. This is in agreement with previous studies (e.g., Padoan et al., 1997; Passot & Vázquez-Semadeni, 1998; Federrath et al., 2008; Price et al., 2011; Burkhart & Lazarian, 2012; Molina et al., 2012) that found that $x$ is lognormally distributed (and $\ln x$ is Gaussian), with

[TABLE]

In these expressions $\sigma_{x}$ and $\sigma_{\ln x}$ are the standard deviations of the $x$ and $\ln x$ distributions, and where $x$ has a unit mean (by definition) and the mean of $\ln x$ is $\mu_{\ln x}=-(1/2)\sigma_{\ln x}^{2}$. The proportionality constant $b$ depends on the nature of the turbulent driving, and ranges from $1/3$ to $1$ for for pure solenoidal or compressive driving respectively (Nordlund & Padoan, 1999; Federrath et al., 2008, 2010). The solid curves in Fig. 2 are Gaussians with $\sigma_{\ln x}$ as given by Eq. (18) with $b=1/3$, appropriate for our solenoidly driven simulations. The agreement is not perfect due to small deviations from the phenomenological $\sigma-\pazocal{M}_{s}$ relation, given by Eq. (18).

4.2. Line-of-sight-averaged Densities

We now discuss an important distribution that will play a crucial role in determining the HI column density PDF. For a column of length $\ell$ we define the average density along a line-of-sight (LOS)

[TABLE]

where $0\leq\ell\leq L_{\rm box}$, and where $L_{\rm box}$ is the simulation box length. For any given $\ell$, different sightlines have different density profiles, and thus the set $x_{\ell}$ form a random variable. We refer to the distribution of $x_{\ell}$ as the “LOS averaged density distribution”.

For turbulent cascade the density is correlated over all scales up to the driving scale. However, the correlation decreases with increasing spatial separation (Vazquez‐Semadeni & Garcia, 2001, hereafter VG01). To obtain an analytic description for the $x_{\ell}$ distribution, we assume that the correlation may be described with a single parameter, $L_{\rm dec}$, hereafter the “decorrelation scale”, such that for ${\ell}<L_{\rm dec}$ the density is effectively constant (i.e., maximally correlated) and for $\ell\geq L_{\rm dec}$ the density cells are uncorrelated. The number of independent density cells along a LOS of length $\ell$ is then

[TABLE]

and the $x_{\ell}$ distribution may be viewed as the sampling distribution of the mean, for which

[TABLE]

This distribution is often encountered in the calculation of errors in repeated measurements (Barlow, 1989). For $\ell\ll L_{\rm dec}$, the LOS contains a single fluctuation $\pazocal{N}\approx 1$, and $\sigma_{x_{\ell}}=\sigma_{x}$. For ${\ell}\gg L_{\rm dec}$, $\pazocal{N}\gg 1$, the LOS contains many turbulent fluctuations, and $\sigma_{x_{\ell}}\ll\sigma_{x}$ as the fluctuations are averaged out.

Fig. 3 shows the ratio of the standard deviations, $\sigma_{x_{\ell}}/\sigma_{x}$, as a function of ${\ell}/L_{\rm box}$, as calculated for our $\pazocal{M}_{s}=2$ and 4.5 simulations. We consider sightlines along the $X$ (solid), $Y$ (dashed), and $Z$ (dotted) directions. The solid blue curve is a fit for the predicted relation

[TABLE]

with the best fitted parameter

[TABLE]

(equivalent to 40 out of 512 cells). Evidently, our simplified treatment for the density correlation gives a reasonable estimate for the $x_{l}$ dispersion.

In Fig. 4 we show PDFs of $\ln x_{\ell}$ (shaded), for $\pazocal{M}_{s}=2$ and 4.5 and ${\ell}/L_{\rm box}=0.25$, 0.5 and 1. The PDFs have distorted Gaussian shapes, becoming narrower with increasing ${\ell}$, as expected from Eq. (22). The blue curves are Gaussians with standard deviations according to Eqs. (17) and (20-23). We conclude that $x_{\ell}$ is indeed well described by a lognormal with $\sigma_{x_{\ell}}=b\pazocal{M}_{s}/\sqrt{\pazocal{N}}$. We note that this relationship derived for $\sigma_{x_{\ell}}$ is complimentary to the column density variance - $\pazocal{M}_{s}$ relationship derived in Burkhart & Lazarian (2012) (their Equation 4 with $A=0.11$) in the limit that ${\ell}=L_{\rm box}$. However the relationship present here is more general and provides a method to determine the driving scale and decorrelation scale via measurement of column density variance.

It is instructive to write $L_{\rm dec}$ in terms of the driving scale. For all our simulations, $L_{\rm drive}=L_{\rm box}/2.5$, and with Eq. (23) we get

[TABLE]

The decorrelation scale is smaller than, but of order of the driving scale. This is because the driving process introduces density (and velocity) correlations, which cascade down to smaller scales. Eq. (24) is a general relation for $L_{\rm dec}$ and $L_{\rm drive}$, although the prefactor may depend on the driving details (e.g. compressional versus solenoidal). VG01 and Fischera & Dopita (2004) also studied the $L_{\rm dec}-L_{\rm drive}$ relation (using alternative methods) and obtained $L_{\rm dec}=0.33L_{\rm drive}$ and $L_{\rm dec}=0.13L_{\rm drive}$, respectively. Our value for the $L_{\rm dec}-L_{\rm drive}$ relation is also in good agreement with that of Kowal et al. (2007).

While $L_{\rm dec}<L_{\rm drive}$, it is typically larger than the sonic scale. For example, Eqs. (11) and (24) suggest that $L_{\rm dec}\geq L_{s}$ as long as $\pazocal{M}_{s}\geq 2.2$. This is important for our model, because $L_{\rm dec}$ represents the scale below which the density becomes effectively uniform. But if $L_{\rm dec}<L_{s}$, then $L_{\rm dec}$ should be everywhere replaced with $L_{s}$.

5. ${\rm HI}-to-{\rm H_{2}}$ in Turbulent Gas

In this section we present atomic and molecular density profiles and integrated HI column density distributions for non-homogeneous, turbulent gas. We use the density field obtained from our MHD simulations, and assume a unidirectional UV flux incident of the box from one side. As we discuss in the Appendix, our results depend weakly on the geometry of the radiation field (e.g. beamed versus isotropic). We solve Eq. (1) to obtain the atomic and molecular density profiles along each LOS. We then integrate the HI densities to obtain the HI column density PDF.

5.1. Basic Parameters

For homogeneous gas, two parameters fully characterize the ${\rm HI}/{\rm H_{2}}$ equilibrium problem, (i) the $\alpha G$ parameter which is proportional to $I_{\rm UV}/n$, and (ii) the dust absorption cross section $\sigma_{g}$, or equivalently $\phi_{g}Z^{\prime}$. Since $n$ is no longer a constant when turbulent fluctuations are present, we define

[TABLE]

where we have replaced $n$ with the volume average $\langle n\rangle$ in Eq. (4). We consider a wide range of $\overline{{\alpha G}}$ values, from the weak ($\overline{{\alpha G}}\ll 2$) to the strong ($\overline{{\alpha G}}\gg 2$) field limits. For $\sigma_{g}$ we assume the standard value $\sigma_{g}=1.9\times 10^{-21}\ {\rm cm^{2}}$ corresponding to $\phi_{g}Z^{\prime}=1$.

The HI column is accumulated over a typical length of

[TABLE]

For a turbulent medium, the density fluctuations have typical lengths of the decorrelation scale $L_{\rm dec}$ (§4.2). Thus, for turbulent medium the ratio

[TABLE]

enters as an additional parameter. The $L_{\rm dec}$-to-$L_{{\rm HI}}$ ratio has the physical meaning of a mean dust opacity over the decorrelation width, denoted by $\tau_{\rm dec}$. The ratio $L_{\rm dec}/L_{{\rm HI}}=\tau_{\rm dec}$ further determines the characteristic number of fluctuations along the HI length, through

[TABLE]

(i.e. Eq. 20 with ${\ell}=L_{{\rm HI}}$), which then controls the ${\rm HI}/{\rm H_{2}}$ structure. Following Eqs. (24,27), the driving scale is related to $\tau_{\rm dec}$ through

[TABLE]

5.2. Profiles

We scale our simulations such that the average optical depth over the box is $\tau_{\rm box}\equiv\sigma_{g}\langle n\rangle L_{\rm box}=5.7$ ensuring HI-to-H2 conversion for all sightlines (since $\tau_{\rm box}\gg 1$). Following Eq. (23) and (28), $\tau_{\rm box}$ sets $\tau_{\rm dec}$ and $\pazocal{N}(L_{{\rm HI}})$, giving $\tau_{\rm dec}=0.45$, $\pazocal{N}(L_{{\rm HI}})=3.2$. We use the above scaling for our results in this section and in §5.3. Following Eqs. (13) and (24), $\tau_{\rm dec}=0.45$ implies $t_{\rm chem}/t_{\rm turb}=0.44(\phi_{g}/\phi_{R})\pazocal{M}_{s}$. Thus, for the highest Mach number we consider ($\pazocal{M}_{s}=4.5$) the chemical time may exceed the turbulent time, unless the H2 formation efficiency is enhanced ($\phi_{R}>2$), or the dust absorption efficiency is reduced ($\phi_{g}<1/2$).

In the upper panels of Fig. 5 we show the density profiles, $x\equiv n/\langle n\rangle$, for three arbitrary sightlines for the $\pazocal{M}_{s}=0.5$, 2 and 4.5. The cloud depth in the abscissa is represented by the mean opacity $\langle\tau({\ell})\rangle\equiv\sigma_{g}\langle n\rangle{\ell}$, ranging from 0 to $\tau_{\rm box}=5.7$. The horizontal bars in each panel represent the decorrelation opacity width $\tau_{\rm dec}=0.45$, which is comparable to the typical scales of density fluctuations. The middle and lower panels show the calculated HI profiles $x_{{\rm HI}}\equiv n_{{\rm HI}}/n$, and the integrated HI column densities, $N_{{\rm HI}}({\ell})\equiv\int_{0}^{l}n_{{\rm HI}}d{\ell}^{\prime}$, for the corresponding LOS, assuming $\overline{{\alpha G}}=2.0$ and $\phi_{g}Z^{\prime}=1$. For comparison, the dashed curves show the exponential decay of $x_{{\rm HI}}$ and the gradual buildup of $N_{{\rm HI}}({\ell})$ as obtained by the uniform-density solution.111In our notation, $N_{{\rm HI}}({\ell})$ refers to the integrated HI column density (that depends on cloud depth), whereas $N_{{\rm HI}}\equiv\lim_{l\rightarrow\infty}N_{{\rm HI}}(l)$ denotes the total (asymptotic) HI column density.

For the subsonic simulation ($\pazocal{M}_{s}=0.5$), the density remains nearly homogeneous and the ${\rm HI}/{\rm H_{2}}$ profiles and the integrated HI column densities remain close to the homogeneous density solution. As the Mach number increases and exceeds unity, density fluctuations become substantial, and the HI (and ${\rm H_{2}}$) density profiles become highly distorted. For the highly supersonic case ($\pazocal{M}_{s}=4.5$) the HI profiles exhibit an extreme scatter, differing by orders of magnitude in some locations. For example, at cloud depth $\langle\tau({\ell})\rangle=3$, $x_{{\rm HI}}\approx 10^{-5}$, $2\times 10^{-2}$, and $0.8$ for the red, yellow and blue LOS. This increasing scatter with $\pazocal{M}_{s}$ is further reflected in $N_{{\rm HI}}({\ell})$ and in the total (asymptotic) HI column density, $N_{{\rm HI}}$. However, because $N_{{\rm HI}}$ is an integrated quantity, the perturbations are (partially) averaged out and the scatter is much smaller than for $x_{{\rm HI}}$. For example, for the three LOS of $\pazocal{M}_{s}=4.5$, the scatter in $N_{{\rm HI}}$ is less than 0.4 dex. In §5.3 we show the calculated PDFs for the total HI columns.

5.3. The HI column density distribution

We integrate the HI profiles along all of the lines-of-sight for the $\pazocal{M}_{s}=0.5$, 2 and 4.5 simulations, and obtain the HI column density distributions. In the upper panels of Fig. 6 we show PDFs of $\log_{10}N_{{\rm HI}}$ for $\overline{{\alpha G}}=0.2$, 2, and 20. The lower panels show the median (solid curve), mean (dotted curve, which almost converges with the median), and the 68.3, 95.5, 99.7 percentiles (about the median) as functions of $\overline{{\alpha G}}$. Evidently, the distributions become wide with increasing $\pazocal{M}_{s}$, and with decreasing $\overline{{\alpha G}}$. For example, for the $\pazocal{M}_{s}=4.5$ simulation the width of the 68.3 percentile is 0.2 dex and 0.5 dex for $\overline{{\alpha G}}=10$ and 0.1 respectively. For the $\pazocal{M}_{s}=0.5$ simulation the width of the 68.3 percentile ranges from 0.03 to 0.1 dex for $\overline{{\alpha G}}=10$ to 0.1. The widening of the HI PDF with increasing $\pazocal{M}_{s}$ reflects the increasing spread of the density PDF, as illustrated in Fig. 2 and in the $\sigma-\pazocal{M}_{s}$ relation, Eq. (17). The HI PDFs become narrow at large $\alpha G$ because in the strong field limit the HI-to-H2 is very sharp as the radiation is absorbed by dust (exponential attenuation), resulting in a weak dependence of the HI column on gas density (Eq. 2).

The median and mean are above the homogeneous solution for small $\overline{{\alpha G}}$ and below it for large $\overline{{\alpha G}}$. Importantly, for all $\overline{{\alpha G}}$, the deviations from the homogeneous solution remain small. For example, for $\overline{{\alpha G}}$ ranging from 0.1 to 10, the deviation is 0.08 to 0.24 dex. for the $\pazocal{M}_{s}=4.5$ case, and is 0.02 to 0.13 dex for $\pazocal{M}_{s}=2$. Thus, the mean HI column is well approximated by the S14 formula for homogeneous gas, as given by Eq. (2), with $n=\langle n\rangle$.

Interestingly, the shapes of the HI PDFs deviate from a lognormal. For all supersonic simulations, the $N_{{\rm HI}}$ distributions are strongly truncated at the highest ends, and have extended tails at the lower-end of the distribution. This is due to the interaction of the propagating radiation with the non-uniform gas and the effect on the H2 self-shielding. At any point inside the cloud, the H2 and HI fractions depend on the local volume density of the gas. However, the H2 self-shielding introduces a dependence on the accumulated H2 column, from the edge to the point of interest. This introduces a non-linear dependence of the ${\rm HI}$/H2 fractions on the volume density profile of the gas. Any positive density perturbation along the column results in a disproportional increase of the H2 fraction and reduced HI fraction, from that point onward. This effect introduces a bias towards lower values of $N_{{\rm HI}}$.

5.4. $\sigma_{\ln N_{{\rm HI}}}$ versus $\overline{{\alpha G}}$, $\pazocal{M}_{s}$ and $\tau_{\rm dec}$

The width of the $\ln N_{{\rm HI}}$ distribution depends on three parameters, $\overline{{\alpha G}}$, the Mach number $\pazocal{M}_{s}$, and the decorrelation opacity $\tau_{\rm dec}$. The first encapsulates H2 formation versus destruction, the second determines the width of the density distribution, and the last determines the frequency of density fluctuations along an HI column. We vary $\tau_{\rm dec}$ by modifying the simulation scaling $\tau_{\rm box}$ (see Eq. 23), and consider $\tau_{\rm dec}$ from 0.01 to 8.

In Fig. 7 we plot $\sigma_{\ln N_{{\rm HI}}}$ (natural logarithm) in the $\overline{{\alpha G}}-\tau_{\rm dec}$ plane, for the $\pazocal{M}_{s}=0.5$, 2, and 4.5 simulations. The dashed line marks the $\tau_{\rm dec}=0.45$ value used for Figs. 5-6. As discussed in §5.3, $\sigma_{\ln N_{{\rm HI}}}$ increases with increasing $\pazocal{M}_{s}$ or with decreasing $\overline{{\alpha G}}$. However, Fig. 7 shows that $\sigma_{\ln N_{{\rm HI}}}$ also has a strong dependence on $\tau_{\rm dec}$. For small $\tau_{\rm dec}$, the number of density fluctuations along a LOS, $\pazocal{N}$, is large ($\pazocal{N}=1+\tau_{\rm dec}^{-1}$). Each LOS then samples a large portion of the (same) parent density PDF, and the different sightlines become more alike. As $\tau_{\rm dec}$ increases, $\pazocal{N}$ decreases, until finally for $\tau_{\rm dec}\gg 1$, $\pazocal{N}\rightarrow 1$. Each LOS is then correlated over the entire $L_{{\rm HI}}$ scale, and its density is approximately uniform, drawn from the parent density PDF. The width of the HI column density is then maximal and reflects the width of the parent density PDF.

In the following section we derive analytic formula for $\sigma_{\ln N_{{\rm HI}}}$ as a function of $\overline{{\alpha G}}$, $\pazocal{M}_{s}$ and $\tau_{\rm dec}$.

6. Analytic Approximation

As is shown in Fig. 5, each LOS has a unique density profile, with a complicated HI density structure. To obtain a simple analytic representation for the distribution of HI columns, we approximate each LOS as containing a uniform density equal to the average along the HI length. I.e. for each LOS we set the density equal to $\langle n\rangle x_{L_{{\rm HI}}}$ (i.e. Eq. 19 with ${\ell}=L_{{\rm HI}}$). Hereafter we use a shortened notation and omit the subscript HI in $x_{L_{{\rm HI}}}$. For each LOS, the HI column density is then given by

[TABLE]

The PDF of $\ln N_{{\rm HI}}$ is thus

[TABLE]

where $\mathrm{d}f/\mathrm{d}\ln x_{L}$ is the PDF of $\ln x_{L}$. Following the discussion in §4.2, $x_{L}$ and $\ln x_{L}$ are approximately Gaussian and lognormal with standard deviations

[TABLE]

respectively. where the second equality follows from Eqs. (17) and (28).

Fig. 8 shows the $\ln N_{{\rm HI}}$ PDFs as obtained from the simulations (also shown in Fig. 8), along with the analytic PDFs, as given by Eq. (31). The locations and the widths of the analytic PDFs roughly follow the PDFs from the simulations, with some differences. First, the shapes of the analytic PDFs are symmetric whereas the simulated PDFs are truncated at the high end and have left-tails. This difference is expected because the analytic approximation does not account for the interaction of the radiation with the density perturbations along each LOS, and as discussed in §5.3, this interaction introduces a preference for small HI columns.

The lower panels of Fig. 8 show the median (solid), mean (dotted), and the 68.3, 95.5, 99.7 percentiles about the median (shaded regions), as functions of $\overline{{\alpha G}}$. The dashed curves are the homogeneous solutions, for comparison. The trend of an increasing dispersion with increasing $\pazocal{M}_{s}$ or with decreasing $\overline{{\alpha G}}$ is in agreement with the numerical results shown in Fig. 6. The median and mean HI columns remain close to the homogeneous solution, also in agreement with the numerical results. An analytic expression for the median $N_{{\rm HI}}$ is obtained by substituting the median

[TABLE]

into Eq. (30), giving

[TABLE]

For moderate $b\pazocal{M}_{s}$ (or if $\tau_{\rm dec}\ll 1$), $x_{L,{\rm med}}\approx 1$ and the median $N_{{\rm HI}}$ remain close to the homogeneous solution (see also Fig. 8). Then, to a good approximation

[TABLE]

similar to Eq. (2) for homogeneous gas, but with ${\alpha G}$ replaced by $\overline{{\alpha G}}$. This result is also confirmed by our numerical computations shown in Fig. 6.

The standard deviation of $\ln N_{{\rm HI}}$ may be approximated by

[TABLE]

where the approximation becomes increasingly accurate for small $\sigma_{\ln N_{{\rm HI}}}$ values. Plugging in Eq. (30) and (33) we get

[TABLE]

In this expression, the nominator is $\sigma_{\ln x_{L}}$ (Eq. 33), introducing the dependence on the turbulence parameters, $b$, $\pazocal{M}_{s}$ and $\tau_{\rm dec}$. As expected, $\sigma_{\ln N_{{\rm HI}}}$ increases with increasing $b\pazocal{M}_{s}$. $\sigma_{\ln N_{{\rm HI}}}$ increases with $\tau_{\rm dec}$, and becomes independent of $\tau_{\rm dec}$ once $\tau_{\rm dec}\gg 1$. For $\tau_{\rm dec}\rightarrow 0$, $\sigma_{\ln N_{{\rm HI}}}\propto\tau_{\rm dec}\rightarrow 0$. The dependence on the radiation intensity enters through the $\overline{{\alpha G}}$ parameter (Eq. 25). For $\overline{{\alpha G}}\ll 1$ (the weak field limit), $\sigma_{\ln N_{{\rm HI}}}\simeq\sigma_{\ln x_{L}}$ and the dispersion is maximal and is independent of $\overline{{\alpha G}}$. For $\overline{{\alpha G}}\gg 1$ (the strong field limit), $\sigma_{\ln N_{{\rm HI}}}\simeq\sigma_{\ln x_{L}}/\ln(\overline{{\alpha G}}/2)$, and the distribution becomes narrow with increasing $\overline{{\alpha G}}$.

Fig. 9 shows $\sigma_{\ln N_{{\rm HI}}}$ as a function of $\tau_{\rm dec}$ and $\overline{{\alpha G}}$, as given by Eq. (38). Like the numerical results (shown in Fig. 7), the standard deviation increases as (a) $\pazocal{M}_{s}$ increases, (b) as $\tau_{\rm dec}$ increases, and (c) as $\overline{{\alpha G}}$ decreases. Deviations from the numerical results exist, and are expected given that the analytic model introduces simplifying assumptions, (a) the density correlations are described by a single decorrelation scale, $\tau_{\rm dec}$, (b) the density distribution is lognormal and follows the $\sigma_{x}-\pazocal{M}_{s}$ relation (Eq. 17), and (c) the ansatz that for each LOS the HI column is given by Eq. (30). The advantage of the analytic approximation is that it provides a smooth solution for $\sigma_{\ln N_{{\rm HI}}}$ as a function of $\pazocal{M}_{s}$ and $b$.

7. Applications to Observations

In this section we present a brief example demonstrating how our results for the width and mean of the HI column PDF may be used to analyze 21 cm observations toward molecular clouds, setting constrains on the Mach number and turbulence driving scale.

Based on 21 cm emission lines from the GALFA-${\rm HI}$ Survey (Peek et al., 2011), Lee et al. (2012) obtained an HI map for the Perseus molecular cloud. Burkhart et al. (2015) derived the HI PDF and found that it to be very narrow, with $\sigma_{N_{{\rm HI}}}/\langle N_{{\rm HI}}\rangle=0.13$. Based on absorption line data from Stanimirović et al. (2014), Burkhart et al. (2015) obtained the Mach number distribution for the cold neutral medium (CNM) around Perseus. They find that $\pazocal{M}_{s}$ ranges from $\pazocal{M}_{s}=1$ to 11, with a median value $\pazocal{M}_{s}=4$.

Unlike the Mach number, the HI PDF, being observed in emission, contains contributions from both the WNM and the CNM phases. The former being typically subsonic, and the latter supersonic (Heiles & Troland, 2003; Wolfire et al., 2003). We decompose the observed $\sigma_{N_{{\rm HI}}/\langle N_{{\rm HI}}\rangle}=0.13$ into CNM and WNM components,

[TABLE]

where the subscripts $C$ and $W$ refer to CNM and WNM, and $\phi_{C}$ and $\phi_{W}$ are the gas mass fractions in these phases. Stanimirović et al. (2014) obtained that around Perseus, $\phi_{C}$ ranges between 0.1 and 0.5, with a median $\phi_{C}=0.35$. Assuming the median $\phi_{C}=0.35$, and assuming that the WNM is subsonic and thus has a negligible HI dispersion (see Figs. 6-9 above), we obtain $\sigma_{N_{{\rm HI},C/\langle N_{{\rm HI},C}\rangle}}=0.22$, or equivalently222For a lognormal distribution, $\sigma_{\ln N_{{\rm HI}}}=\ln^{1/2}(1+\sigma_{N_{{\rm HI}}/\langle N_{{\rm HI}}\rangle}^{2})$. $\sigma_{\ln N_{{\rm HI},C}}=0.22$, for the CNM in Perseus.

Inspecting the $\pazocal{M}_{s}=4.5$ panel in Fig. 7, we see that for $\overline{{\alpha G}}$ within the CNM range $(\overline{{\alpha G}})_{\rm CNM}=1-8$ (see §2), the $0.22$ contour is obtained for $\tau_{\rm dec}=L_{\rm dec}/L_{{\rm HI}}=0.06$ to $0.3$, respectively. These values of $\tau_{\rm dec}$ correspond to $L_{\rm drive}/L_{{\rm HI}}=0.3-1.5$ (Eq. 29), i.e. the driving scale is of order of the HI scale-length. Assuming typical CNM density, $\langle n\rangle\approx 30$ cm*-3*, and for $\phi_{g}Z^{\prime}=1-2$ as suggested by Lee et al. (2012) for Perseus, we obtain (with Eq. 26) $L_{\rm drive}\sim 1-8$ pc.

If we do not neglect the WNM contribution in Eq. (39), $\sigma_{\ln N_{{\rm HI},C}}$ would be smaller, further reducing $\tau_{\rm dec}$ and $L_{\rm drive}$. Our numerical computations in Fig. 7 are based on the MHD simulations that assume solenoidal driving ($b=1/3$). For compressional driving or mixed driving the contours in Fig. 7 would shift downwards (since increasing $b$ is similar to increasing $\pazocal{M}_{s}$; Eq. 17) and $\tau_{\rm dec}$ and $L_{\rm drive}$ will again decrease. Thus, our derived driving scale is an upper limit.

A similar conclusion may be drawn from the analytic model. Eq. (38) predicts the standard deviation of $\ln N_{{\rm HI}}$ as a function of $b$, $\pazocal{M}_{s}$, $\tau_{\rm dec}$ and $\overline{{\alpha G}}$. Plugging in $\sigma_{\ln N_{{\rm HI},C}}=0.22$ and $\overline{{\alpha G}}=1-8$, and inverting Eq. (38) we get

[TABLE]

where $C=0.28-0.47$ for $\overline{{\alpha G}}=1-8$ respectively. This relation between the Mach number, the driving forcing parameter $b$ and the driving scale ($\tau_{\rm dec}=0.2L_{\rm drive}/L_{{\rm HI}}$; Eq. 29) is shown in Fig. 10 for $b=1/3$, 1/2 and 1, corresponding to solenoidal, mixed and compressional drivings. The width of the strips correspond to the width of the $\overline{{\alpha G}}=1-8$ CNM range. For the median $\pazocal{M}_{s}=4$ (dashed line), $L_{\rm drive}/L_{{\rm HI}}$ ranges within 0.2-0.7, 0.1-0.3, and 0.02-0.07 for $b=1/3$, $b=1/2$ and $b=1$.

Interestingly, narrow HI PDFs were recently reported for more Galactic clouds: Ophiuchus, Orion A, Orion B, California, MonR2 and Rosette (Imara & Burkhart, 2016), suggesting that the presence of a small driving scale, of order of, or smaller than the HI length, might be a general feature in molecular clouds and/or a considerable WNM component to the HI column density.

However, there are caveats to this analysis. First, the decomposition of the HI PDF from emission line measurements into CNM and WNM components, which requires CNM/WNM fraction from absorption measurements, is uncertain. Thus, it would be valuable to use absorption line measurements directly to infer the HI PDF of the CNM. Second, accurate measurements of the Mach number are difficult due to line of sight blending of features in position-position-velocity space. Third, the width of the HI PDF may be also affected by (a) variations in $I_{\rm UV}$ across the observed region, (b) the chosen cutoff criteria of the observed HI cube; the velocity range, and the spatial extent, (c) the finite angular resolution which may smooth out density fluctuations, and (d) optically thick ${\rm HI}$ that may produce artificially narrow HI PDFs (Burkhart et al., 2013b). While for Perseus optical thickness is probably not a major issue (Lee et al. 2015 obtain $\sim 10\%$ corrections for the optically thick gas), it might be important in other GMCs. A thorough observational analysis that addresses these uncertainties will be provided elsewhere. Model limitations are discussed in §8.

8. Summary and Conclusions

In this paper we have studied the HI-to-H2 transition, and the ${\rm HI}$ column densities, maintained by far-UV radiation, in turbulent media. We have used a suite of MHD simulations to produce realistic turbulent density distributions (§4). The density field is nearly lognormal and the dispersion follows the PDF variance-sonic Mach number ($\sigma-\pazocal{M}_{s}$) relation (§4). We find that for supersonic gas, the density decorrelation length is related to the driving scale through $L_{\rm dec}=0.2L_{\rm drive}$. The decorrelation length determines the number of density fluctuations along a LOS, through $\pazocal{N}=1+{\ell}/L_{\rm dec}$, where $\ell$ is the length-scale along the LOS. This simple description, allows us to model the scale dependent LOS averaged density distribution (§4.2) and the HI column density PDF analytically. We note that while pure MHD simulations are scale-free, when the depth dependent HI-to-H2 transition is of interest (or any other chemical species), the adopted simulation box length plays an important role. The adopted box scale determines $L_{\rm dec}$, which in turn defines the number of density fluctuations along an HI column ($\pazocal{N}=1+L_{{\rm HI}}/L_{\rm dec}$), thus affecting the ${\rm HI}/{\rm H_{2}}$ structure (see §5.1 and §5.4).

As we show in §5.2, once the turbulence becomes supersonic, strong density fluctuations are developed, and the atomic-to-molecular density profiles are significantly distorted relative to those for homogeneous uniform density media. As a result, different lines of sight (LOS) differ in their total accumulated HI column density, $N_{{\rm HI}}$. We calculate the probability density function (PDF) for $N_{{\rm HI}}$ (§5.3-5.4), as a function of the governing physical parameters: (a) the sonic Mach number $\pazocal{M}_{s}$, (b) the effective dissociation parameter $\overline{{\alpha G}}$ and (c) the decorrelation opacity $\tau_{\rm dec}$, which is related to the turbulence driving scale ($\tau_{\rm dec}\equiv L_{\rm dec}/L_{{\rm HI}}=0.2L_{\rm drive}/L_{{\rm HI}}$).

We find that the mean and median $N_{{\rm HI}}$ are affected by turbulence, but as long as $b\pazocal{M}_{s}\sim 1$, or if $\tau_{\rm dec}\ll 1$, $N_{{\rm HI}}$ may be well approximated by the S14 and BS16 formula for uniform-density gas,

[TABLE]

where

[TABLE]

is the effective dimensionless dissociation parameter, where $I_{\rm UV}$ is the free-space intensity of the far-UV field and $\langle n\rangle$ is the volume density. Here $\sigma_{g}$ is (cm2) is the dust-grain absorption cross section per hydrogen nuclei averaged over the Lyman-Werner dissociation band.

The major effect of turbulent density fluctuations, is in producing a spread in the HI column distribution. For subsonic gas the density is nearly uniform and the HI PDF is very narrow and the solutions converge to the S14 and BS16 formula for uniform-density gas. As the Mach number increases, density fluctuations becomes substantial and the HI PDF also widens. As discussed in §5, the HI PDF also depends on $\overline{{\alpha G}}$ and $\tau_{\rm dec}$, becoming wider for small $\overline{{\alpha G}}$ or small $\tau_{\rm dec}$. In §6, we present an analytic formula, Eq. (38), for the standard deviation of the HI PDF as a function $\pazocal{M}_{s}$, $\overline{{\alpha G}}$, and $\tau_{\rm dec}$.

We demonstrate how our model may be combined with 21 cm observations toward GMCs to constrain turbulent parameters (§7). For Perseus (and in other Galactic clouds), the very narrow observed HI PDF may suggest small-scale driving, potentially pointing to the importance of multi-scale turbulent driving in the CNM (Haverkorn et al., 2008; Yoo & Cho, 2014). Alternatively, the narrow PDF may be caused by small-scale decorrelation lengths induced by the abrupt change in the chemical and thermal properties at the HI-to-H2 transition. Observational caveats are discussed in §7.

Our results are for single-phased gas, irradiated by far-UV radiation. In §7 we describe how our results may be still applied to a mixed CNM/WNM medium, given estimates of the CNM and WNM fractions. This marks the importance of observations of both emission and absorption 21 cm lines, which provide constrains on the CNM fraction and its Mach number.

The values of observed HI column densities may depend on geometry (i.e. slabs versus spheres, or slab inclination) and the number of HI-to-H2 transition layers along the observed sightlines. However, importantly, the turbulent parameters, $\pazocal{M}_{s}$ and $\tau_{\rm dec}$ are obtained from the standard deviation of the logarithmic HI column, $\sigma_{\ln N_{{\rm HI}}}$ (Eq. 38). Thus, our analysis is robust against any multiplication of the HI column, including, number of clouds, geometry factors, and inclination corrections.

A basic assumption made in our analysis is that the HI and H2 are in chemical steady state (see §3 for a discussion of the timescales). The coupling of time dependent chemistry and turbulence may affect the HI PDF in various ways. For example, rapidly changing density fluctuations may alter the dispersion in the HI distribution, since HI and H2 in different LOS do not have time to react to the fast density changes. Turbulent mixing, may increase the mean HI column by transferring molecular gas from inner shielded to outer unshielded regions, where it can rapidly dissociate and form ${\rm HI}$. We plan to investigate time-dependent effects in a future work.

We conclude that the atomic-to-molecular density profiles and the HI column density are affected by the turbulent nature of the CNM. For moderate Mach numbers, the homogeneous solution still provides a good estimate for the mean HI column density. The standard deviation of the HI PDF contains a useful information regarding the turbulence properties. Our model, combined with 21 cm observations, may be used to constrain the sonic Mach number and turbulence driving scale of cold atomic gas.

We thank Enrique Vazquez-Semadeni, Sahar Shahaf and Ewine van Dishoeck for fruitful discussions. We thank the referee for constructive comments that improved our paper. This work was supported in part by the PBC Israel Science Foundation I-CORE Program grant 1829/12, by DFG/DIP grant STE1869/2-1 GE625/17-1, and by the Raymond and Beverly Sackler Tel Aviv University - Harvard/ITC Astronomy Program. B.B. acknowledges support from the NASA Einstein Postdoctoral Fellowship and the Raymond and Beverly Sackler TAU-ITC Visiting Researcher fund.

Appendix A Beamed versus isotropic irradiation

In our radiative transfer analysis we have assumed unidirectional beamed irradiation. Beamed irradiation is expected in the proximity of strong far-UV sources. For clouds immersed in the diffuse Galactic interstellar radiation field, isotropic irradiation may be a better approximation. S14 showed that for uniform-density clouds, exposed to isotropic fields the HI column density is given by

[TABLE]

where $\langle\mu\rangle\equiv 0.8$ is a geometrical factor. For a given ${\alpha G}$, the flux (normal to the slab) for isotropic radiation is half the beamed flux, and hence the HI column density is smaller (see §2.3 in S14 for a full discussion and derivation of Eq. A1). Comparing Eq. (A1) with the beamed solution (Eq. 2), we see that

[TABLE]

where $N_{{\rm HI}}^{\rm iso}$ and $N_{{\rm HI}}^{\rm bm}$ are the HI columns produced by isotropic and beamed fields, respectively. Eq. (A2) provides a simple transformation between the beamed and the isotropic solutions, for uniform-density gas. Such a transformation is particularly useful given that an isotropic field calculation is much more time consuming.

However, Eq. (A2) may be less accurate for turbulent gas, because the density fluctuations will have different effects on the ${\rm HI}/{\rm H_{2}}$ structures, for the different field geometries. As an illustrative example, consider a diffuse medium, and a very dense clump. Let the clump be a thin disk of radius $R$, located near the slab surface at $z=r=0$ (cylindrical coordinates), where $z$ is the axis normal to the slab along which we integrate $N_{{\rm HI}}$. For beamed radiation, all rays originating at $r\leq R$ will pass through the clump and $N_{{\rm HI}}$ would be very small for these sightlines, whereas for LOS with $r>R$, the gas is diffuse and $N_{{\rm HI}}$ would be very large. For isotropic radiation, inclined rays (originating at $r>R$) would penetrate behind the clump and thus $N_{{\rm HI}}(r\leq R)$ would increase compared to the beamed case. On the other hand, for the isotropic field, gas at $r>R$ may be still partially shielded by the clump leading to a decrease in $N_{{\rm HI}}(r>R)$ compared to the beamed case. These two opposite effects will tend to smooth out the contrast in the $N_{{\rm HI}}$ map, and also reduce the width of the HI PDF.

To test the effects of isotropic versus beamed irradiation, we have carried out a calculation of the ${\rm HI}/{\rm H_{2}}$ structure for isotropic radiation impinging the $\pazocal{M}_{s}=4.5$ simulation box, with $\overline{{\alpha G}}=2.6$ (the mean CNM value). The resulting HI column density map is shown in Fig. 11 (left panel). The corresponding PDF is shown in the right panel (red shaded). For comparison, the middle panel, and the gray histogram in the right panel, show $0.8\times N_{{\rm HI}}$ for beamed irradiation with $\overline{{\alpha G}}=2.6/1.6=1.63$ (as suggested by the scaling of A2). We see that the two maps and PDFs are not identical, as would be the case for a uniform-density medium. For the isotropic case, the HI map is smoother and the low-$N_{{\rm HI}}$ features are more extended, as expected. The width of the isotropic PDF is somewhat narrower ($\sigma_{\ln N_{{\rm HI}}}=0.15$), and has a lower mean ($\langle N_{{\rm HI}}\rangle=1.5\times 10^{20}$ cm*-2*) compared to the beamed field ($\sigma_{\ln N_{{\rm HI}}}=0.18$, $N_{{\rm HI}}=2.1\times 10^{20}$ cm*-2*).

Importantly, while there are differences between the beamed and isotropic cases, especially when considering individual lines of sight, the relative differences in the mean and standard deviations are small ($28\%$ and $18\%$, respectively). We conclude that to a good approximation, Eq. (A2) may be used as a transformation from our (efficiently computed) beamed field results to isotropic radiation fields.