Chemical Abundances in a Turbulent Medium -- H$_2$, OH$^+$, H$_2$O$^+$, ArH$^+$
Shmuel Bialy David Neufeld Mark Wolfire Amiel Sternberg Blakesley Burkhart





1Harvard-Smithsonian Center for Astrophysics, 60 Garden street, Cambridge, MA, USA
2Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 20910, USA
3University of Maryland, Astronomy Department, College Park, MD 20742, USA
4School of Physics & Astronomy, Tel Aviv University, Ramat Aviv 69978, Israel
5Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave., New York, NY, 10010, USA
6Max-Planck-Institut für extraterrestrische Physik (MPE), Giessenbachstr., 85748 Garching, FRG, Germany
7Department of Physics and Astronomy, Rutgers University,  136 Frelinghuysen Rd, Piscataway, NJ 08854, USA
^⋆sbialy@cfa.harvard.edu
ArH^+, OH^+, H_2O^+ in a Turbulent Medium
Bialy et al.
The Astrophysical Journal. Submitted: Jul. 18, 2019.  Accepted: Sept. 25, 2019


Supersonic turbulence results in strong density fluctuations in the interstellar medium (ISM), which have a profound effect on the chemical structure. 
Particularly useful probes of the diffuse ISM are the ArH^+, OH^+, H_2O^+ molecular ions, 
which are highly sensitive to fluctuations in the density and the H_2 abundance. 
We use isothermal magnetohydrodynamic (MHD) simulations of various sonic Mach numbers, ℳ_s, and density decorrelation scales, y_ dec, to model the turbulent density field.
We post-process the simulations with chemical models and obtain the probability density functions (PDFs) for the H_2, ArH^+, OH^+ and H_2O^+ abundances.
We find that the PDF dispersions  increases 
with increasing ℳ_s and y_ dec, as the magnitude of the density fluctuations increases, and as they become more coherent.
Turbulence also affects the median abundances:
when ℳ_s and y_ dec are high, low density regions with low H_2 abundance become prevalent, resulting in an enhancement of ArH^+ compared to OH^+ and H_2O^+.
We compare our models with Herschel observations. 
The large scatter in the observed abundances, as well as the high observed ArH^+/OH^+ and ArH^+/H_2O^+ ratios are 
naturally reproduced
 by our supersonic (ℳ_s=4.5), large decorrelation scale (y_ dec=0.8) model, supporting a scenario of a large-scale turbulence driving.
The abundances also depend on the UV intensity, CR
ionization rate and the cloud column density, and the observed scatter may be influenced by fluctuations in these parameters.



















§ INTRODUCTION








A diverse collection of molecules and molecular ions has been detected in the interstellar medium (ISM).
The abundances of observed species and their ratios are often used to constrain cloud properties: the temperature, incident ultraviolet (UV) radiation flux, cosmic-ray ionization rate (CRIR), and the density of H nuclei, n ≡ n( H) + 2n( H_2).

Diffuse clouds in the ISM are commonly analyzed by chemical models for photodissociation-regions (PDRs; e.g., ).
<cit.> used constant-n  PDR models to study the abundances of ArH^+, OH^+, and H_2O^+, and compared their model predictions with observations.
These molecular-ions are sensitive to the H_2 abundance, which in turn depends on the 
 the UV intensity, the CRIR, the cloud visual extinction, and the gas density.



 showed that for all sightlines, there is a discrepancy between their theoretical models and the observations, 
as the models always under-predict the ArH^+ abundance relative to OH^+ and H_2O^+. 

To explain this discrepancy, they were forced to invoke a two-cloud population model in which low A_ V clouds were the source for the ArH^+, while higher A_ V clouds produced the OH^+ and H_2O^+.
Furthermore, the observations show scatter in the column density ratios N( ArH^+)/N( H), N( OH^+)/N( H), and N( H_2O^+)/N( H), for different lines-of-sight (LoS).
 attributed the scatter to variations in the CRIR.


A common simplifying assumption in PDR models is that the gas density in the cloud is constant (or is smoothly varying with cloud depth).
However, this assumption may be wrong.
Observations suggest that molecular clouds as well as cold atomic clouds exhibit turbulent supersonic motions, as is evident from, e.g.: the observed superthermal
linewidths, the line-width size relation, the fractal structures of molecular clouds, and studies of the power spectrum/bispectrum <cit.>. 






The turbulence, when supersonic, produces strong density fluctuations in the gas, which alter the chemical structure. 

<cit.> have shown that in a turbulent medium, density fluctuations 
perturb the H and H_2 abundances in the gas,
resulting in large dispersions in the atomic columns, N( H),
for different LoS. 
The variance in the N( H) distribution was related to the governing turbulent parameters: the sonic Mach number, , and the decorrelation scale of the density field, L_ dec, which is proportional to the turbulence driving scale, L_ drive.


In this paper we extend upon the  and   analyses. 
We study how turbulence affects the abundances of H_2, ArH^+, OH^+, and H_2O^+, and compare our models with observations.
We find that our turbulent-cloud models may explain the observed scatter in the observed abundances, as well as the high ArH^+/OH^+ and ArH^+/H_2O^+ ratios.
We demonstrate how the observations may be used to constrain the turbulent parameters,  and L_ drive.











The paper is organized as follows.
In  <ref>  we review the chemical network and discuss the governing physical parameters.

In <ref> we discuss the model ingredients, including our MHD simulations and the chemical PDR models. 
In <ref> we present our results for the abundance distributions in turbulent clouds.
In <ref> we compare our models
with observations.
We discuss the limitations of the model and our conclusions in <ref> and <ref>.



























§ CHEMISTRY





The formation and destruction pathways include two-body reactions, surface chemistry on dust-grains, and cosmic-ray and UV photoreactions.

We define the abundance of species i, x( i) ≡ n( i)/n, where n(i) is the density of species i, and n is the density of hydrogen nuclei.




 §.§ H and H_2



The H_2 is formed on the surfaces of dust grains.
The removal of H_2 is through photodissociation by Lyman-Werner (LW) radiation (11.2-13.6 eV),
and by cosmic-ray (CR) ionization.
The H and H_2 steady-state abundances
obey

    x( H_2)/x( H) =  Rn/D_0 f_ att+ζ_ tot ,






and by element conservation, x( H) + 2x( H_2) ≃ 1.
Here, R= 3 × 10^-17 cm^3 s^-1

is the
 H_2 formation rate coefficient for 
 formation on dust-grains
(the value is for a 100 K gas),
and D_0=5.8 × 10^-11 I_ UV s^-1 
<cit.>, 

where
I_ UV is the interstellar radiation field in units of the value given by <cit.>.
The attenuation factor, f_ att, accounts for the absorption of LW photons by dust and in H_2 lines.
For normal incidence radiation, penetrating a slab on two sides 



    f_ att = 1/2f_ sh(N^l( H_2)) e^-σ_g N^l  + 1/2f_ sh(N^r( H_2)) e^-σ_g N^r .

Here σ_g =1.9 × 10^-21 cm^2  is the dust absorption cross section per hydrogen nucleus <cit.>, and
f_ sh is the H_2 self-shielding function 
<cit.> 
for which we adopt b=2 km/s (see <ref> and appendix <ref>).
N and N( H_2) are the column densities of hydrogen nuclei and H_2, from cloud edge to the point of interest, and the superscripts l and r denote integration from the left and right cloud edges, respectively.





In Eq. (<ref>), ζ_ tot≈ 2.3 × 1.9 ×ζ is the total CR removal rate of H_2 and where ζ is the primary CRIR of H.
The factor 2.3 translates from primary CR ionization per H to total (primary+secondary) ionization per H_2 
<cit.>, and the additional factor 1.9 accounts for H_2 destruction by chemical reactions with molecular ions which are produced by CR ionization (i.e., the H_2 abstraction reactions below)[



We have fixed the 2.3 and 1.9  factors at values appropriate for cloud interiors.
The value 1.9 is found by taking an average value 

from our modeling results accounting for the H_2 destruction from all chemical pathways.

Towards cloud edge, these factors vary, however
in these regions H_2 destruction is in any case dominated by UV photodissociation, not by CR processes.
].
We define ζ_-16≡ζ/(10^-16 s^-1)






















 §.§ ArH^+, OH^+, and H_2O^+



Here we discuss the basic formation and destruction chemistry
for ArH^+, OH^+, and H_2O^+ (see also
 and ).
The production of ArH^+ is initiated by CR ionization followed by hydrogen abstraction 

    Ar +  cr   → Ar^+ +  e , 
     Ar^+ +  H_2   → ArH^+ +  H .

The formation sequence is moderated by the recombination of Ar^+ with polycyclic aromatic hydrocarbons (PAH), and PAH^- <cit.>.
Dissociative recombination of ArH+ is unusually slow
(), and its destruction is dominated by
proton transfer
with H_2

    ArH^+ +  H_2→ Ar +  H_3^+ .

If the H_2 abundance is low, x( H_2)≲ 10^-4, the removal of ArH^+ proceeds mainly through charge transfer with O and C.




The production of OH^+ 
is initiated by CR ionization of H,  followed
by charge transfer and H abstraction,

    H +  cr   → H^+ +  e ,     
     H^+ +  O   → O^+ +  H ,  
     O^+ +  H_2   → OH^+ +  H .

If the H_2 abundance is sufficiently high, OH^+ is formed through the sequence

    H_2+ cr   → H_2^+ +  e , 
     H_2^++ H_2   → H_3^+ +  H , 
     H_3^++ O   → OH^+ +  H_2 .

Further abstraction reactions with H_2 destroy OH^+ and lead to the formation of H_2O^+ and H_3O^+

    OH^+ +  H_2   → H_2O^+ +  H , 
     H_2O^+ +  H_2   → H_3O^+ +  H .

OH+ and H2O+ are also destroyed  by dissociative recombination and photodissociation (see Fig. 2 in
).

Because of the abstraction reactions with H_2, the OH^+ and H_2O^+ abundances are very sensitive to x( H_2), down to very low H_2 abundances, x( H_2) ∼ 10^-3 - 10^-4. 
At still
lower H_2 abundances, OH^+ and H_2O^+ are formed via solid-state chemistry on dust grains
<cit.>, and their abundances are independent of x( H_2). 
In this
regime, atomic hydrogen on grains reacts with atomic oxygen to form OH. The OH
can be desorbed from the grain surface or react with another hydrogen to form 
H_2O. 
Upon desorption, OH and H_2O undergo a charge 
transfer   
with H^+ forming OH^+ and H_2O^+.












In Fig. <ref> we show the ArH^+, OH^+ and H_2O^+ abundances as functions of x( H_2), as computed by our PDR models (see <ref> below)
for I_ UV=1, ζ_-16= 4, and for various densities, n.
For a fixed x( H_2) the abundances of all the molecular ions increase with decreasing n. 
The formation rate is proportional to ζ, while destruction is proportional to n, so all three molecular ions have abundances that increase with ζ/n.

At a constant n, the x( ArH^+) is a non-monotonic function of x( H_2): it increases with x( H_2) when x( H_2) is very small, peaks at x( H_2) ∼ 10^-4, and then decreases with x( H_2).
This behavior may be understood as follows. 
When x( H_2) is very small, the abstraction reaction (<ref>) is slow, and is heavily interrupted by 
Ar^+ recombinations with PAHs.
As a result, the formation rate is ∝ζ multiplied by a branching ratio ∝ x( H_2).
Destruction, on the other hand, is dominated by charge transfer with O and C, and is independent of x( H_2). Thus the abundance scales with x( H_2).
At high x( H_2), the abstraction reaction proceeds rapidly such that every Ar ionization leads to the formation of ArH^+, and the formation is independent of x( H_2).
However, removal is now dominated by interactions with H_2.
Thus the overall abundance scales inversely with x( H_2). 
Similar considerations explain the non-monotonic behavior of the OH^+ and H_2O^+ abundances.





§ MODEL INGREDIENTS







 §.§ Chemical PDR Models



We ran a large series of chemical PDR models similar to those described by <cit.> and <cit.>.
Each model is characterized by (1) I_ UV/n, (2) ζ/n, and (3) total visual extinction, A_ V, through the cloud (denoted A_ V (tot) by ).


The models span the parameter space I_ UV/n=20 to 10^-3 cm^3, ζ/n= 10^-20 to 4 × 10^-16 cm^3 s^-1, and A_ V=3× 10^-4 to 8 mag.
For any given set of parameters, the model computes the temperature and molecular abundances (assuming steady-state) as functions of cloud depth, hereafter the z coordinate.

Generally the fractional densities inside the cloud, x(i), depend on the four parameters, I_ UV/n, ζ/n, the total cloud visual extinction A_ V, and cloud depth z.
However, for ArH^+, OH^+, and H_2O^+, the dependence on z and A_ V is fully captured by the variation of the H_2 fraction, x( H_2).
To this end, 
we re-express the molecular-ion abundances as functions of x( H_2), reducing the rank of the parameter-space by 1  ({z,A_ V}→ x( H_2)). 
We construct lookup tables for x(i) (i= ArH^+, OH^+, H_2O^+) as functions of I_ UV/n, ζ/n and x( H_2).




















 §.§ MHD Simulations










To model the density field in the turbulent medium, we solve the ideal MHD equations

    ∂ρ/∂ t + ∇· (ρ)    = 0  , 
    ∂ρ/∂ t + ∇·[ ρ + ( p + B^2/8 π)  I - 1/4 π B B]    =  f ,  
    ∂ B/∂ t - ∇× (× B)    = 0  .

Here ρ is density,

B is magnetic field, p is the gas pressure, I is the identity matrix and f is the specific force.
We use a third-order-accurate hybrid essentially nonoscillatory scheme on a 3D grid of 512^3  cells.
The B field is initially aligned along the y axis.
We assume 
 periodic boundary conditions, and an isothermal equation of state p = c_s^2 ρ.
By constrained transport we numerically construct ∇· B = 0.
For the source term f, we assume a random large-scale solenoidal driving at a wave number k≈ 2.5 (i.e. 1/2.5 the box size). 
For more details see <cit.>.



Each simulation is defined by two parameters, the sonic and Alfvénic Mach numbers, ≡⟨ || ⟩/c_s, ℳ_A ≡⟨ || ⟩/ _A, where  is the velocity, c_s  and _A are the isothermal sound speed and the Alfvén speed, and ⟨·⟩ denotes averages over the entire simulation box.
We consider three simulations, =0.5, 2.0 and 4.5, corresponding to subsonic, transonic, and supersonic medium. All simulations have ℳ_A = 2.
As discussed in , ℳ_s strongly affects the variance of the density field,
 and hence the H-H_2 structure.
 However, the H-H_2 structure is only weakly sensitive to the Alfvénic number.
As the simulations do not include gravity and chemistry the density and size are scale-free 
<cit.>.
We adopt = 6.3 pc for the simulation box length, and ⟨ n ⟩ = 30 cm^-3, corresponding to a mean cloud column ⟨ N ⟩≡⟨ n ⟩ = 5.8 × 10^20 cm^-2 or ⟨ A_ V⟩=0.3 mag.



In Figure <ref> we show density slices, and density PDFs, for the =0.5 and 4.5 simulations.
The density is nearly homogeneous in the subsonic simulation, while in the supersonic case shocks develop and lead to strong density fluctuations in the gas.



The n-PDF is close to a lognormal, 
with a dispersion that increases with the sonic Mach number,
as σ_ln (n/⟨ n ⟩)^2 ≃ln [1+ (b )^2], where
b=1/3 is the forcing parameter <cit.>.











The large-scale driving yields a density field that is correlated over large scales.
The decorrelation scale
 in the simulation is ≡ L_ dec/ = 0.08, i.e., ≈ 40 cells (). 
L_ dec is the characteristic scale over which
correlations in the density field decrease:
For averaging over length L, the density variance decreases with increasing L, and 
 L_ dec is the characteristic scale for this decreasing trend
 (see 4.2 in . cf.  for alternative definitions).
In supersonic gas, density fluctuations are still present on scales below L_ dec, all the way down to the sonic-length, L_s (see footnote 10  below). 
As evident in Fig. <ref>, L_ dec is indeed a typical length-scale for the 
large-scale density modes
 in the simulations.
The value of  has a strong influence on the resulting chemical structure (see <ref> below).

To study the effects of variations in , we crop each of the simulations after 51 cells, producing thin slabs of thickness 1/10 of the original boxes.
The relative decorrelation length in the slabs (along the short axes) is ≃ 0.8. The slab thickness is still set to
= 6.3 pc, as for the original simulation boxes.








 §.§ Numerical Procedure














As our focus is on the effects of turbulence, we fix the non-turbulent parameters
I_ UV=1, ζ_-16=4, ⟨ n ⟩  = 30 cm^-3 and ⟨ A_ V⟩  = 0.3, and study the dependence on the turbulent parameters  and .
The effects of I_ UV, ζ and A_ V are studied in . 
The LoS direction from the observer to the cloud (simulation) along which we integrate the column is aligned with the z axis.
z is also the axis along which the radiation propagates (i.e., as we assume slab geometry).
In the Appendix we explored other LoS orientations, along x or y, and find that the results are weakly sensitive to the LoS direction.



We post process each of the simulations as follows:


  * The 3D density field, n(x,y,z), is used as an input to the chemical model.

  * For each LoS, (x,y), we solve Eqs. (<ref>-<ref>)[This requires iterations since at any point x( H_2) depends on the entire H_2 profile on both sides (through the H_2 self-shielding function).
Clouds of sufficiently large column density may be modeled by a combination of two semi-infinite slabs, for which the H and H_2 profiles may be obtained analytically by integrating over the H_2 shielding function <cit.>.
However, this approximation is justified only if N( H) ≫ 2/σ_g ln[D_0G/(2Rn)+1] ≈ 5 × 10^20 cm^-2 
where G ≈ 3× 10^-5 <cit.>.
Our clouds do not satisfy this condition and thus we must solve for the two-sided irradiation.] on
a fine logarithmic grid, and  then interpolate the results back onto the simulation box.

  * Given x( H_2)(x,y,z) and n(x,y,z), we use the lookup tables to 
obtain the ArH^+, OH^+ and H_2O^+ abundances in all the cells, x(i)(x,y,z).

  * We also integrate the abundances
along the LoS direction, z, yielding the species column densities for all LoS, N(i)(x,y).

We now have the abundance and column density PDFs for the six cases:
=(0.5, 2, 4.5) for =(0.08, 0.8).




		
		
		













§ RESULTS















In this section we present the probability distribution functions (PDFs) of the abundances x(i), and the columns, N(i).






















 §.§ Abundance PDFs



To demonstrate the effect of the density fluctuations on the chemical structure, in Fig. <ref> we plot 
the joint-PDFs for n, x( ArH^+), x( OH^+) and x( H_2O^+) with x( H_2), and their marginal PDFs, for the subsonic =0.5 (left) and the supersonic 4.5 (right) simulations (both with =0.08).
In the subsonic case the density PDF is narrow, with most cells having density n≈⟨ n ⟩ = 30 cm^-3.
The H_2 distribution exhibits a spread, from log x( H_2) ∼ -5 to -1,  corresponding to regions near cloud boundaries and center.
The ArH^+, OH^+ and H_2O^+ PDFs as functions of x( H_2) follow the shapes of the constant density n=30 cm^-3 models (see Fig. <ref>), as they should.

In the supersonic case, strong density fluctuations are developed, and the n PDF is wide.
The onset of H_2 self-shielding leads to a rapid growth in x( H_2) with increasing cloud depth (column density) resulting in a bimodal PDF.
At a given depth, higher gas densities result in more efficient H_2 formation and thus a positive correlation between x( H_2) and n.
At cloud boundaries, the correlation is linear, with x( H_2) = Rn/(D_0/2).
This is evident as the thin diagonal density enhancement extending from log x( H_2), log n = (-3, 3) to (-6, 0).
For n ≲ 10 cm^-3, there is an additional density enhancement, which also extends diagonally, but at an offset of factor of 2 to the left.
This corresponds to LoS that are so diffuse that radiation from the far side of the cloud penetrates, resulting in a factor of 2 higher dissociation rate.  

The density fluctuations in the =4.5 simulation enhance the dynamic range of x( H_2) and thus the dynamic ranges of ArH^+, OH^+ and H_2O^+.
At a given x( H_2), the density fluctuations also lead to a spread in the ArH^+, OH^+ and H_2O^+ PDFs, through their dependencies on ζ/n.





Since x( H_2) and n are correlated, the PDF peaks as functions of x( H_2) no longer follow any of the constant density models.
For example, as x( H_2) increases from 10^-4 to 0.5, the ArH^+ falls from ≈ 10^-9 to ≲ 10^-13, which correspond to the n=30 and 
300 cm^-3
models in Fig. <ref>, respectively.
  
The marginal PDFs of the molecular ions are composed of a broad and a narrow component.
In the case of OH^+ and H_2O^+ the narrow components correspond to formation on dust near cloud boundaries (where x( H_2) is small) which lead to characteristic OH^+ and H_2O^+ abundances, that depend weakly on x( H_2) or n.
In the case of ArH^+, the narrow PDF reflects the maximum in the x( ArH^+) versus x( H_2) relation (which occurs at x( ArH^+) ∼ 10^-9 for ⟨ n ⟩ =30 cm^-3).




 §.§ Column Density PDFs




 In Fig. <ref>
we show the column density  (or column density ratios) PDFs


in the form of violin charts (where the width is proportional to the PDF dispersion), 
as obtained 
for all of our simulations, =0.5, 2 and 4.5, for =0.08 (left) and 0.8 (right).

All PDF widths increase with increasing 
as a result of the density fluctuations.
The widths vary among species as is determined by the individual dependencies of x(i) on n (see <ref>).
The column PDFs are further affected by the density as cells with higher density 
contribute more to the column, N(i) = ∫ x(i) n d z.

For example, the N( H) PDFs are all very narrow because x( H) is anti-correlated with n, thus the cells where x( H) is high have a low contribution to N( H), and vice versa.
By contrast, x( H_2) correlates with n, and thus the N( H_2) PDFs are very wide.










The medians also show a dependence on .
The median ArH^+ increases with , both at =0.8 and 0.08.
This is because the density field is characterized by a lognormal PDF in which low densities are more common than high (i.e., n_ med<⟨ n ⟩).
These diffuse regions are deficient in H_2, and thus enhanced in ArH^+, leading to an overall increase in the median ArH^+ value.
This is evident by comparing the =0.5 and 4.5 abundance PDFs in Fig. <ref>.
For =0.8, the median OH^+ and H_2O^+ decrease with increasing , as the H_2-poor gas 
depresses OH^+ and H_2O^+ formation.

At small , the trends with  are less pronounced as there are many density fluctuations along the sightlines, which, upon summation, partially cancel out (see <ref> below).
The sum is not evenly weighted as cells with
high density contribute more to the total integrated column. 
This is the reason why the median OH^+ trend is reversed at =0.08:
While x( OH^+) is depressed in the low x( H_2) (low n) regions, the depression is sufficiently weak so that the x( OH^+)-rich (high n) regions are able to (over)compensate through the larger weight they have in the integrated column.























  §.§.§ The Role of 
 



The PDFs in Fig. <ref>  depend strongly on the decorrelation scale,  
(see discussion in <ref>).
Most prominent is the increase in the PDF dispersions as  increases.
This is because when  is large, the density fluctuations are coherent along the LoS (the number of fluctuations ∼ 1/).
Since the density still varies between LoS, there are large variations in the abundances among the sightlines and the PDFs are wide.
On the other hand, when  is small, 
each LoS contains several (∼ 10) density fluctuations along it.
Upon integration, positive and negative fluctuations average-out, resulting in narrower PDFs.
The  H_2 PDF is particularly wide in the =0.8 case, because 
at any point in the cloud, the H_2 is affected by all cells along the LoS through H_2 self-shielding and thus coherent density fluctuations strongly affect H_2.

Generally, because at large  the density field is more coherent, the column PDFs more closely resemble the x(i) PDFs, in this limit.
For example, in the =0.8 case 
the OH^+ and H_2O^+ PDFs both have peaks at low values (at log N( i)/N( H) ≈ -9.3 and -11.1), which are absent in the =0.08 case.
These peaks correspond to low-density clouds with low x( H_2),
in which the OH^+ and H_2O^+ are formed through the dust-catalysis formation routes, reflecting the narrow peaks in the x(i) PDFs in Fig. <ref>.
For =0.08, these features disappear as the low-x(i) regions are averaged over with  regions with higher x(i), where the high-x(i) regions have a larger weight.




The dependence of PDF width on  may be understood analytically for the case of the N PDF (also in the case of N( H), see ).
Following , 
we approximate the density field in two cells along the LoS as (a) two uncorrelated variables if Δ z>L_ dec where Δ z is the seperation along the LoS, or (b) assume full correlation if Δ z < L_ dec .
The number of fluctuations along a LoS is then

    𝒩 = 1+ 1/ .

In this approximation, the column density is the sum of 𝒩 independent variables, and its standard deviation obeys

    σ_N/⟨ N ⟩≃σ_n/ ⟨ n ⟩/√(𝒩)≃b /√(1+ 1/) .

The PDF dispersion increases with  and with increasing fluctuation size,   (or with decreasing 𝒩).





 Eq. (<ref>) is in 10-30% agreement with our numerically calculated standard deviations for all six simulations.




In the limit ≫ 1,
σ_N/⟨ N ⟩→σ_n/⟨ n ⟩ independent of .
In this limit, each cloud corresponds to a single density fluctuation in the large-scale density field.
The column PDF of species may then be obtained directly from the uniform  
density models, weighting each model with the n-PDF.





The decorrelation scale is related to the driving scale:

    L_ dec =  ϕ L_ drive ,

or in dimensionless units, 
= ϕ (L_ drive/L_ cloud).
 In our simulations ϕ=0.2
 (see 4.2 in ),  in agreement with analytic and numerical studies which find ϕ≈ 0.1-0.3[
 For a line-width size relation with exponent 1/2, the sonic scale L_s = L_ drive/^2 <cit.>, and is   typically ≪ L_ dec. For example, for =4.5, L_s/L_ drive=0.05, while L_ dec/L_ drive=0.2.
 ] 
 <cit.>.
  
The dependence of the species PDFs on  may be potentially used as a method to constrain the turbulence driving scale in clouds (). 
















§ COMPARISON TO OBSERVATIONS








 §.§ Observations


The three molecular ions discussed here, OH^+, H_2O^+ and ArH^+, have all been observed extensively in the diffuse Galactic ISM.   Absorption-line observations of their submillimeter rotational transitions near 972, 1115 and 618 GHz, respectively, have been carried out by the Herschel Space Observatory towards several regions of massive star formation that serve as background continuum sources.  Such observations typically reveal multiple absorption components arising in diffuse foreground gas that is unassociated and spatially-separated from the continuum sources; thanks to the differential rotation of the Galaxy, multiple diffuse clouds along each sight-line may be distinguished kinematically.  At the typical density in the diffuse ISM, these molecular ions are found primarily in the ground rotational state; thus, the observations yield robust estimates of the molecular column densities that do not depend on precise knowledge of the gas density or of the rate coefficients for collisional excitation.   Because OH^+ and H_2O^+ have rotational lines that show hyperfine structure, a deconvolution must be performed prior to the determination of their column densities <cit.>. 

From Herschel observations reported <cit.> in previous literature, <cit.> identified fifteen velocity intervals for which all three molecular ions were detected within the spectra of 4 background continuum sources: W31C, G34.3, W49N, and W51e.  The choice of velocity intervals is somewhat arbitrary, in that the observed absorption cannot be reliably decomposed into the contributions from individual clouds; thus each velocity interval may contain multiple diffuse clouds.  The molecular column densities for these velocity intervals, along with ancillary estimates of the HI column densities derived from  21 cm observations <cit.>, comprise the observational data with which we compare the model predictions. 




 §.§ The Grand PDFs


Since the observations are probing through crowded regions in the Galactic plane, 
the sightlines likely contain several clouds superimposed along the LoS <cit.>.
This is also supported by the observed spectra which show several distinct peaks, and is expected 
based on the large 
observed column densities N( H)_ obs=0.84-6.1 × 10^21 cm^-2. 
For our turbulent cloud models ⟨ N( H) ⟩≈ 4 × 10^20 cm^-3
giving m=2-15 clouds along the LoS. 



To compare our model with observations, we autoconvolve our single-cloud column PDFs (<ref>) as follows.

    
  * For each simulation we randomly draw 10^4 LoS, and calculate the various column densities.
    
  * We repeat step (1) m times, and for each LoS we co-add the respective columns.
    This produces the “m-cloud PDF", for m  superimposed clouds along the LoS.
    We do that for any m ∈ [1,2,...,24][Here m_ max=24 is chosen to satisfy m_ max≫ N( H)_ obs, max/⟨ N( H) ⟩≈ 15.].
    
  * We stack all the m-cloud PDFs with equal weighting to produce 
    a single “grand" PDF (for each species), which accounts for the superposition of several (unknown number) of clouds along the LoS.
    
    
    
    


We show the grand PDFs as functions of N( H) in Fig. <ref>
for (a) =0.5, =0.08, (b) =4.5, =0.08, and (c)
=4.5, =0.8.
In each panel, the three strips enclose the 68, 95 and 99.7 percentiles about the median (at constant N( H)).
The observations are shown as dots with errorbars.
Here,

N( H) on the y-axis encodes information regarding the number-of-clouds, with m ≈ 2 N( H)/(10^21  cm^-2).
This is thanks to the fact that the N( H)-PDFs of a single cloud are very narrow (see Fig. <ref>).
On the other hand, from an observational point, N( H) is a direct observable.
On the x-axis, since the columns are normalized to N( H), the median values are independent of m (or N( H)), and all align up with the corresponding single-cloud values.

The dispersions of the PDFs (across the x axis) decrease with increasing m, as 1/√(m), 
(see discussion in <ref>).
The PDF shapes approach Gaussians with increasing m, as required by the central limit theorem.
Importantly, Fig. <ref> allows the comparison of observations with theory, even though the different observational sightlines may have different number of clouds, and thus a different predicted PDF width.


Comparing the three cases in Fig. <ref> we observe the familiar trends as for the single-cloud PDFs (<ref>):
(1) the increase in the PDF dispersions as the gas becomes supersonic and density fluctuations develop,
(2) the further increase in the dispersions as the decorrelation scale increases,  and
(3) the shift of the ArH^+ median to higher values, and the OH^+ and H_2O^+ medians to lower values,  as  and  increase.


The model that explains best the observations is the =4.5, =0.8 model



in which the PDF dispersions are largest, and where the N( ArH^+)/N( OH^+) ratios are highest.
The large decorelation scale combined with Eq. (<ref>) suggests that the density fluctuations are driven on scales L_ drive≫ L_ cloud.
However, even in this model the observational scatter is too large requiring either still larger  and , or additional fluctuations in the non-turbulent parameters as we discuss below.




§ DISCUSSION






    


 §.§ Variations in non-turbulent parameters


    As our focus in this paper is on the effect of the turbulence induced density perturbations, we fixed the “non-turbulent parameters", I_ UV, ζ, ⟨ n ⟩ and ⟨ A_ V⟩, and studied the behavior as function of the “turbulent parameters",  and .
    However, in the real ISM, the non-turbulent parameters may vary between clouds inducing additional variations to the abundances.
    In fact, even in our high- high- model, in which the dispersion is maximal, it is still not large enough to account for the scatter in the observations. 
    This suggests additional fluctuations in the non-turbulent parameters, and/or a higher Mach number and/or higher decorrelation lengths.


    It is important to keep in mind that the abundances do not depend on I_ UV, ζ and ⟨ n ⟩ separately, but instead are determined by the two ratios I_ UV/⟨ n ⟩ and ζ/⟨ n ⟩.
    Therefore, if there exist positive correlations between the UV intensity, CRIR, and density, any fluctuation in each of these parameters will result in only a small fluctuation in their ratios.  
    Indeed, under thermal steady-state conditions, 
    the density of the cold neutral medium (CNM, T ∼ 100 K)  is predicted to be positively correlated with I_ UV <cit.>. 
    ζ and I_ UV may also correlate given that both are produced by 
    massive stars (the latter in the supernovae remnants, at the stars' death), however, this  would depend on the CRs diffusion scale.
    
    Additional independent observational constraints of these parameters, their fluctuations, and correlations, are of great interest 
and will allow to reduce the degeneracy with the turbulent parameters.
    In the absence of additional measurements, the observed dispersion may still be used to place upper limits on the turbulent parameters and on the dispersion in the non-turbulent parameters.
    
    
        
    
    
    
    
    



 §.§ Clouds and Density Fluctuations
   

    What defines “a cloud" in a turbulent ISM?
    Generally, one may define a cloud as a region that exceeds some threshold density, or mass.
    However, with this definition the cloud characteristics: its mass, size and internal structure, would depends on the (somewhat arbitrary) density/mass threshold value.
    Instead, in this study we defined the cloud based on the UV (more particularly, the LW) radiation field.
    We explicitly assumed that each cloud is irradiated by the mean  radiation field, 
    such that each cloud self-shields itself and does not affect the other clouds.
    For this condition to be fulfilled, the clouds need to be sufficiently separated in space (or in velocity space).
    Such a description of the ISM, as an ensemble of externally irradiated clouds, is useful but 
    is clearly an over-simplification of the real ISM.
    
    It is also worth reiterating the importance of taking into account the superposition of clouds along the LoS   (<ref>).
    While the mean and median abundances relative to H are the same for a single cloud and several superimposed clouds, the dispersion is sensitive to the number of clouds, and it decreases as ≃ 1/√(m).
    While for m=1 the shapes of the PDFs are highly irregular, we find that as m increases the PDFs approach Gaussians, as predicted by the central limit theorem.
    


 §.§ Thermal Phases of the ISM
 
  
In this study we used isothermal MHD simulations, which mimic a single phase medium (the CNM).
In the real ISM, cooling/heating processes introduce a thermal instability which results in a 
multiphasic medium:
 the warm and cold neutral media (WNM, CNM), T ≈ (6000,60) K, ⟨ n ⟩≈ (0.3,30) cm^-3 <cit.>, where intermediate temperatures are thermally unstable  <cit.>.
MHD simulations which include heating-cooling processes find a bimodal PDFs of  T and n, with a non-negligible mass within the unstable region <cit.>, where the bimodality is less pronounced in strongly turbulent gas
<cit.>.


How a bimodal PDF may affect the chemistry?
This critically depends on the mixing-scale of CNM/WNM structures.
Let λ_ c be a characteristic  CNM scale, and consider two limiting cases:
(a) λ_ c≪ L_ cloud, and
(b) λ_ c≳ L_ cloud,
where L_ cloud is defined so that it is externally irradiated by LW radiation (see <ref>).
In case (a), the CNM and WNM are mixed on small scales, and each cloud contains both phases, i.e., the entire bi-modal PDF is sampled in each cloud.
In this case, there would be more  H_2-poor gas compared to the pure CNM case, 
which in turn would affect the abundances of other species.
On the other hand, for case (b), the scale is large, so that each cloud samples only a part of the bimodal PDF, e.g., just the higher CNM densities. 
Our study thus resembles case (b). 
This discussion highlights the importance 
of determining the scales of WNM-CNM structures in the ISM
<cit.>.

Even if an individual cloud is described by only the CNM phase,
if L_ drive≫ L_ cloud,   the driving
may take place in a multiphase CNM-WNM medium.
Both the L_ drive-L_ dec and the σ_n/⟨ n ⟩- relations  are derived from idealized isothermal simulations.
In future studies, it is important to generalize these relations to the more realistic case of driving in a multiphase medium.



 §.§ The H_2 self-shielding function
   

The H_2 self-shielding function of <cit.> and <cit.> depend on the 
Doppler broadening parameter, b.
The assumption is that the gas is turbulent on small scales (compared to the length-scale over which the H_2 abundance change), aka, micro-turbulence.
Micro turbulence broadens the H_2 lines, in a similar manner to thermal broadening.
In appendix <ref> we study the dependence on b and show that this dependence is weak.

However, turbulence induces gas motions that may be correlated on various scales, up to the driving scale, and the micro-turbulence assumption may fail.
 <cit.> have suggested an alternative shielding function that approximates this effect, however it still lacks the information on the velocity correlations, and the density fluctuations.
 On the other hand direct computations of H_2 shielding are possible given the velocities and densities in the simulation cells, however, this is computationally expansive as it requires radiative transfer calculations with high spectral resolution.
The effects of a turbulent velocity field on the H_2 column are expected to be important for clouds of small total column (for I_ UV=1 and n=3-300 cm^-3, N ∼ 10^18-10^20 cm^-2, see Fig. 8 in
 <cit.>, and for large velocity dispersions, with gas motions that are coherent over large scales.
 We plan to investigate the question of H_2 line absorption in turbulent gas elsewhere, with the goal of deriving an all-purpose self-shielding function.
 



 §.§ Comparison to Previous Work


In a previous study,  studied the atomic-to-molecular transition and the HI column PDF. 
They presented numerical results and  provided analytic formula for the N( H) dispersion as a function of and  (and the effective dissociation parameter, see their Eq. 38).
For similar reasons to those discussed here (<ref>), 
the dispersion in N( H) increases with  and .
  applied their model  
to 21 cm observations of HI columns towards the Perseus molecular.
As the N( H) values towards Perseus are relatively uniform <cit.>,  concluded that   is small.
Conversely, in this study we find that the dispersions  in the abundances of ArH^+, OH^+, and H_2O^+, over the various LoS are relatively large, requiring a large .

The large dispersion observed here which is not seen in Perseus, may suggest that non-turbulent parameters, i.e., ζ or I_ UV, contribute considerably to the observed scatter in the  abundances.
However, a direct comparison  is complicated as the 
observations in the two studies probe different kind of clouds, and in different Galactic environments.
Furthermore, other factors also affect the dispersion, such as the number of clouds along the LoS, and the multiphase structure of the ISM (see 7-8 in   for a further discussion). 

Our adopted value, ζ_-16 = 4, optimizes the fit of our supersonic high  model to the measured column densities of OH^+, H_2O^+, and ArH^+.
Our value is a factor of (1.4-2) higher than that derived previously by .
For the subset of 15 sources considered here (in  which ArH^+ was detected), and using their two-cloud population model with uniform densities, they derived ⟨logζ/n_50⟩ = -15.25 with a standard error on the mean ± 0.07 (see their Table 2).
For n=30 cm^-3 this corresponds to ⟨ζ_-16⟩≈ (2.0-2.8),
where we included a downward correction to  account for the higher I_ UV/n assumed here, following the relation ζ∝ (I_ UV/n)^0.7 (see 3.4 in ).
For a further comparison to <cit.>, see Table 2 in .




§ CONCLUSIONS


We studied how turbulence-induced density fluctuations affect the chemical structure of diffuse clouds, our conclusions are:


  * ArH^+, OH^+, and H_2O^+ are sensitive probes of the 
H_2 abundance, which is in turn very sensitive to density fluctuations, because of H_2 self-shielding. 

  * When < 1, the density fluctuations are weak, and the resulting abundances converge to those predicted by the former uniform-density models.
For supersonic gas, the density fluctuations become strong, and as a result the H_2, ArH^+, OH^+ and H_2O^+ abundances in the cloud span a large range.

  * The shapes of the x(i)-PDFs are irregular, typically double-peaked, and are determined through the unique chemical formation pathways, e.g., OH^+ and H_2O^+ gas-phase formation versus formation on dust.

  * The column density PDFs become wider as  increases, and also as   increases.
At high , the density field is more coherent, enhancing the effect of the density fluctuations on the chemical structure.
On the other hand, when  is small, positive and negative fluctuations partially average out.

  * The medians are also affected by the density fluctuations.
With increasing  and , the median ArH^+ increases, while OH^+ and H_2O^+ decrease.
When  and  are high, regions with low H_2 abundance are prevalent, which favor ArH^+ over OH^+ and H_2O^+.

  * The observed abundances have a considerable scatter, and high ArH^+-to-OH^+ and ArH^+-to-H_2O^+ ratios.
These suggest supersonic gas with large decorrelation-scales, which in turn correspond to large-driving scale, L_ drive≫ L_ cloud.

  * The abundances also depend on ζ/⟨ n ⟩, I_ UV/⟨ n ⟩ and A_ V, and the observed scatter may result from fluctuations in these parameters.


More independent observations are needed to break the degeneracy between the turbulent and non-turbulent parameters. On the theoretical part, the next steps are to (a) develop more realistic models that include the interplay of turbulence and thermal phases, and (b) a more accurate H_2 shielding function that applies to turbulent medium.
With these advances, observations of molecular ions may be used as a robust tool to constrain the  driving and strength of interstellar turbulence, and/or variations in the interstellar fields: the UV intensity and the cosmic-ray ionization rate.
  


















We thank the anonymous referee for valuable comments and suggestions that improved this work.
We thank the Center for Computational Astrophysics at the Flatiron Institute for hospitality and funding where some of this research was carried out.
This work was also supported by the German Science Foundation via DFG/DIP grant STE 1869/2-1 GE 625/17-1 at Tel-Aviv University.

56
natexlab#1#1

[Arnold et al.(2015)Arnold, Thomas, Loch, Abdel-Naby, &
  Ballance]Arnold2015
Arnold, I., Thomas, E., Loch, S. D., Abdel-Naby, S., & Ballance, C. P. 2015,
  Journal of Physics B: Atomic, Molecular and Optical Physics, 48, 175005

[Audit & Hennebelle(2010)]Audit2010
Audit, E., & Hennebelle, P. 2010, Astronomy and Astrophysics, 76, 1

[Bell et al.(2006)Bell, Roueff, Viti, & Williams]Bell2006
Bell, T. A., Roueff, E., Viti, S., & Williams, D. A. 2006, Monthly Notices of
  the Royal Astronomical Society, 371, 1865

[Bialy et al.(2017a)Bialy, Bihr, Beuther, Henning, &
  Sternberg]Bialy2017b
Bialy, S., Bihr, S., Beuther, H., Henning, T., & Sternberg, A.
  2017a, The Astrophysical Journal, 835, 126

[Bialy et al.(2017b)Bialy, Burkhart, &
  Sternberg]Bialy2017
Bialy, S., Burkhart, B., & Sternberg, A. 2017b, The Astrophysical
  Journal, 843, 92

[Bialy & Sternberg(2015)]Bialy2015a
Bialy, S., & Sternberg, A. 2015, Monthly Notices of the Royal Astronomical
  Society, 450, 4424

[Bialy & Sternberg(2016)]Bialy2016a
—. 2016, The Astrophysical Journal, 822, 83

[Bialy & Sternberg(2019)]Bialy2019
—. 2019, The Astrophysical Journal, 881, 160

[Bialy et al.(2015)Bialy, Sternberg, Lee, Petit, &
  Roueff]Bialy2015c
Bialy, S., Sternberg, A., Lee, M.-Y., Petit, F. L., & Roueff, E. 2015, The
  Astrophysical Journal, 809, 122

[Burkhart et al.(2009)Burkhart, Falceta-Goncalves, Kowal, &
  Lazarian]Burkhart2009
Burkhart, B., Falceta-Goncalves, D., Kowal, G., & Lazarian, A. 2009, The
  Astrophysical Journal, 693, 250

[Chepurnov et al.(2015)Chepurnov, Burkhart, Lazarian, &
  Stanimirovic]Chepurnov2015
Chepurnov, A., Burkhart, B., Lazarian, A., & Stanimirovic, S. 2015, The
  Astrophysical Journal, 810, 33C

[Choi & Stone(2012)]Choi2012
Choi, E., & Stone, J. M. 2012, Astrophysical Journal, 747

[Draine(1978)]Draine1978
Draine, B. T. 1978, The Astrophysical Journal Supplement Series, 36, 595

[Draine & Bertoldi(1996)]Draine1996
Draine, B. T., & Bertoldi, F. 1996, The Astrophysical Journal, 468, 269

[Federman et al.(1979)Federman, Glassgold, & Kwan]Federman1979
Federman, S. R., Glassgold, A. E., & Kwan, J. 1979, The Astrophysical Journal,
  227, 466

[Federrath et al.(2008)Federrath, Klessen, &
  Schmidt]Federrath2008
Federrath, C., Klessen, R. S., & Schmidt, W. 2008, The Astrophysical Journal,
  688, L79

[Field(1965)]Field1965
Field, G. B. 1965, The Astrophysical Journal, 142, 531

[Field et al.(1969)Field, Goldsmith, & Habing]Field1969
Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, The Astrophysical
  Journal, 155, L149

[Fischera & Dopita(2004)]Fischera2004
Fischera, J., & Dopita, M. A. 2004, The Astrophysical Journal, 611, 919

[Gazol & Kim(2013)]Gazol2013
Gazol, A., & Kim, J. 2013, The Astrophysical Journal, Volume 765, Issue 1,
  article id. 49, 8 pp. (2013)., 765

[Glassgold & Langer(1974)]Glassgold1974
Glassgold, A. E., & Langer, W. D. 1974, The Astrophysical Journal, 193, 73

[Gnedin & Draine(2014)]Gnedin2014
Gnedin, N. Y., & Draine, B. T. 2014, The Astrophysical Journal, 795, 37

[Heiles & Troland(2003)]Heiles2003
Heiles, C., & Troland, T. H. 2003, The Astrophysical Journal, 586, 1067

[Heyer et al.(2009)Heyer, Krawczyk, Duval, & Jackson]Heyer2009
Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, The Astrophysical
  Journal, 699, 1092

[Hill et al.(2008)Hill, Benjamin, Kowal, Reynolds, Haffner, &
  Lazarian]Hill2008
Hill, A. S., Benjamin, R. A., Kowal, G., et al. 2008, The Astrophysical
  Journal, 686, 363

[Hollenbach et al.(2012)Hollenbach, Kaufman, Neufeld, Wolfire, &
  Goicoechea]Hollenbach2012
Hollenbach, D., Kaufman, M. J., Neufeld, D., Wolfire, M., & Goicoechea, J. R.
  2012, The Astrophysical Journal, 754, 105

[Imara & Burkhart(2016)]Imara2016
Imara, N., & Burkhart, B. 2016, The Astrophysical Journal, 829, 102

[Indriolo & McCall(2012)]Indriolo2012
Indriolo, N., & McCall, B. J. 2012, The Astrophysical Journal, 745, 91

[Indriolo et al.(2015)Indriolo, Neufeld, Gerin, Schilke, Benz,
  Winkel, Menten, Chambers, Black, Bruderer, Falgarone, Godard, Goicoechea,
  Gupta, Lis, Ossenkopf, Persson, Sonnentrucker, van der Tak, van Dishoeck,
  Wolfire, & Wyrowski]Indriolo2015
Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, The Astrophysical
  Journal, 800, 40

[Kaufman et al.(1999)Kaufman, Wolfire, Hollenbach, &
  Luhman]Kaufman1999
Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, The
  Astrophysical Journal, 527, 795

[Kowal et al.(2007)Kowal, Lazarian, & Beresnyak]Kowal2007
Kowal, G., Lazarian, A., & Beresnyak, A. 2007, The Astrophysical Journal, 658,
  423

[Le Petit et al.(2006)Le Petit, Nehme, Le Bourlot, Roueff,
  et al. Abergel, Abgrall, Abgrall, Bakes, Balakrishnan, Bates, Bell, Black,
  Black, Black, Bluhm, Boger, Bohlin, Burke, Burton, Boissé, Cardelli,
  Casu, Cazaux, Cecchi-Pestellini, Chambaud, Decamp, de Jong, Dickinson,
  Draine, Draine, Dufton, Faure, Federman, Federman, Fitzpatrick, Fitzpatrick,
  Fitzpatrick, Flannery, Flower, Flower, Flower, Flower, Flower, Gredel, Gry,
  Habing, Hollenbach, Jaquet, Jura, Joulain, Kaufman, Kaufman, Kopp, Kwan,
  et al. Lacour, Laor, Launay, Launay, Lavendy, Le Bourlot, Le Bourlot, Le
  Bourlot, Le Bourlot, Le Bourlot, Le Bourlot, Le Bourlot, Le Petit,
  Le Petit, et al. Le Picard, Le Teuff, Lee, Lennon, Mathis, Mathis,
  Mendoza, Meyer, Meyer, Monteiro, Morton, Nehmé, Neufeld, Patriarchi,
  Patriarchi, Péquignot, Péquignot, Péquignot, Prasad, et al.
  Rachford, Roberge, Roueff, Roueff, Roueff, Savage, Schroder, Staemmler,
  Stephens, Sternberg, Sternberg, Teyssier, Tielens, Tielens, et al. Tumlinson,
  Turner, van Dishoeck, Weingartner, Wells, & Wilson]LePetit2006
Le Petit, F., Nehme, C., Le Bourlot, J., et al. 2006, The Astrophysical
  Journal Supplement Series, 164, 506

[Lee et al.(2012)Lee, Stanimirović, Douglas, Knee, Di
  Francesco, Gibson, Begum, Grcevich, Heiles, Korpela, Leroy, Peek, Pingel,
  Putman, & Saul]Lee2012
Lee, M.-Y., Stanimirović, S., Douglas, K. A., et al. 2012, The
  Astrophysical Journal, 748, 75

[McCourt et al.(2018)McCourt, Oh, O'Leary, & Madigan]McCourt2018
McCourt, M., Oh, S. P., O'Leary, R., & Madigan, A. M. 2018, Monthly Notices of
  the Royal Astronomical Society, 473, 5407

[McKee & Ostriker(2007)]McKee2007
McKee, C. F., & Ostriker, E. C. 2007, Annual Review of Astronomy and
  Astrophysics, 45, 565

[Neufeld & Wolfire(2016)]Neufeld2016
Neufeld, D. A., & Wolfire, M. G. 2016, The Astrophysical Journal, 826, 183

[Neufeld & Wolfire(2017)]Neufeld2017b
—. 2017, The Astrophysical Journal, 845, 163

[Piontek & Ostriker(2007)]Piontek2007
Piontek, R. A., & Ostriker, E. C. 2007, The Astrophysical Journal, 663, 183

[Röllig et al.(2006)Röllig, Ossenkopf, Jeyakumar,
  Stutzki, & Sternberg]Rollig2006
Röllig, M., Ossenkopf, V., Jeyakumar, S., Stutzki, J., & Sternberg, A.
  2006, Astronomy and Astrophysics, 451, 917

[Roman-Duval et al.(2010)Roman-Duval, Jackson, Heyer, Rathborne, &
  Simon]Roman-Duval2010
Roman-Duval, J., Jackson, J. M., Heyer, M., Rathborne, J., & Simon, R. 2010,
  The Astrophysical Journal, 723, 492

[Sanchez et al.(2005)Sanchez, Alfaro, & Perez]Sanchez2005
Sanchez, N., Alfaro, E. J., & Perez, E. 2005, The Astrophysical Journal, 625,
  849

[Saury et al.(2014)Saury, Miville-Deschênes, Hennebelle,
  Audit, & Schmidt]Saury2014
Saury, E., Miville-Deschênes, M.-A., Hennebelle, P., Audit, E., &
  Schmidt, W. 2014, Astronomy & Astrophysics, 567, A16

[Schilke et al.(2014)Schilke, Neufeld, Müller, Comito, Bergin,
  Lis, Gerin, Black, Wolfire, Indriolo, Pearson, Meneten, Winkel,
  Sánchez-Monge, Möller, Godard, & Flalgarone]Schilke2014
Schilke, P., Neufeld, A., Müller, H. S. P., et al. 2014, Astronomy &
  Astrophysics, A29

[Sonnentrucker et al.(2015)Sonnentrucker, Wolfire, Neufeld, Flagey,
  Gerin, Goldsmith, Lis, & Monje]Sonnentrucker2015
Sonnentrucker, P., Wolfire, M., Neufeld, D. A., et al. 2015, The
  Astrophysical Journal, 806, 49

[Sternberg & Dalgarno(1995)]Sternberg1995
Sternberg, A., & Dalgarno, A. 1995, The Astrophysical Journal Supplement
  Series, 99, 565

[Sternberg et al.(2014)Sternberg, Petit, Roueff, &
  Bourlot]Sternberg2014
Sternberg, A., Petit, F. L., Roueff, E., & Bourlot, J. L. 2014, The
  Astrophysical Journal Supplement Series, 790, 10S

[Stutzki et al.(1998)Stutzki, Bensch, Heithausen, Ossenkopf, &
  Zielinsky]Stutzki1998
Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998,
  Astronomy and Astrophysics, 336, 697

[Tielens & Hollenbach(2002)]Tielens1985a
Tielens, A. G. G. M., & Hollenbach, D. 2002, The Astrophysical Journal, 291,
  722

[Vazquez-Semadeni(1994)]Vazquez-Semadeni1994a
Vazquez-Semadeni, E. 1994, The Astrophysical Journal, 423, 681

[Vazquez‐Semadeni & Garcia(2001)]VazquezSemadeni2001
Vazquez‐Semadeni, E., & Garcia, N. 2001, The Astrophysical Journal, 557, 727

[Walch et al.(2011)Walch, Wuensch, Burkert, Glover, &
  Whitworth]Walch2011a
Walch, S., Wuensch, R., Burkert, A., Glover, S., & Whitworth, A. 2011, The
  Astrophysical Journal, 733, 47

[Waters & Proga(2019)]Waters2019
Waters, T., & Proga, D. 2019, The Astrophysical Journal, 876, L3

[Winkel et al.(2017)Winkel, Wiesemeyer, Menten, Sato, Brunthaler,
  Wyrowski, Neufeld, Gerin, & Indriolo]Winkel2017
Winkel, B., Wiesemeyer, H., Menten, K. M., et al. 2017, Astronomy &
  Astrophysics, 600, A2

[Wolfire et al.(1995)Wolfire, Hollenbach, McKee, Tielens, &
  Bakes]Wolfire1995
Wolfire, M., Hollenbach, D., McKee, C., Tielens, A., & Bakes, E. 1995, The
  Astrophysical Journal, 443, 152

[Wolfire et al.(2010)Wolfire, Hollenbach, & McKee]Wolfire2010
Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, Astrophysical Journal,
  716, 1191

[Wolfire et al.(2003)Wolfire, McKee, Hollenbach, &
  Tielens]Wolfire2003
Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, The
  Astrophysical Journal, 587, 278







§ A. DEPENDENCE ON THE DOPPLER BROADENING PARAMETER


The Doppler broadening parameter, b, enters the H_2 self-shielding function.
In Fig. <ref> we show the PDF median and dispersion for all the species and for all the simulations considered here, for different b values.
As b increases, the difference in velocities of gas elements along the LoS increases and self-shielding becomes less efficient.
As a result, the H_2 column decreases.
However, this dependence is very weak because absorption by the H_2 line Doppler cores is effective at relatively small H_2 columns N( H_2) ≈ 10^16-10^17 cm^-2, corresponding to total columns N ∼ 10^18-10^20 cm^-2 (for I_ UV=1 and n=3-300 cm^-3; see Fig. 8 in
 <cit.>.
Since the clouds considered here have ⟨ N ⟩ = 5.8 × 10^20 cm^-2, absorption by the Doppler cores is sub-dominant.
Instead, H_2 self shielding is dominated by the Damping wings (i.e., Lorentz broadening) independent of b.
For similar reasons, the H and the total columns are insensitive to b, in agreement with the finding of .
 
The molecular ions, show a stronger dependence on b, especially, ArH^+.
 This is because the ions, and especially ArH^+, 
 are formed efficiently in regions of low x( H_2).
 These regions typically lie closer to the cloud edges where 
Doppler broadening is important. 
In any case, variations in  and  have a stronger effect on the PDFs.


	
		


§ B. DEPENDENCE ON THE LOS ORIENTATION


The magnetic field in the simulations introduces a preferred direction (the  B field is initiated along the y axis).
In Fig. <ref> we plot the PDF median and dispersion for all the species and for all the simulations considered here, for LoS parallel to x, y or z.
As evident from the figure, the results are weakly sensitive to the direction of the LoS.





