HI in Molecular Clouds: Irradiation by FUV plus Cosmic Rays

We extend the analytic theory presented by Sternberg et al. (2014) and Bialy&Sternberg (2016) for the production of atomic hydrogen (HI) via FUV photodissociation at the boundaries of dense interstellar molecular (H$_2$) clouds, to also include the effects of penetrating (low-energy) cosmic-rays for the growth of the total HI column densities. We compute the steady-state abundances of the HI and H$_2$ in one-dimensional gas slabs in which the FUV photodissociation rates are reduced by depth-dependent H$_2$ self-shielding and dust absorption, and in which the cosmic-ray ionization rates are either constant or reduced by transport effects. The solutions for the HI and H$_2$ density profiles and the integrated HI columns, depend primarily on the ratios $I_{\rm UV}/Rn$ and $\zeta/Rn$, where $I_{\rm UV}$ is the intensity of the photodissociating FUV field, $\zeta$ is the H$_2$ cosmic-ray ionization rate, $n$ is the hydrogen gas density, and $R$ is the dust-surface H$_2$ formation rate coefficient. We present computations for a wide range of FUV field strengths, cosmic-ray ionization rates, and dust-to-gas ratios. We develop analytic expressions for the growth of the HI column densities. For Galactic giant molecular clouds (GMCs) with multiphased (warm/cold) HI envelopes, the interior cosmic-ray zones will dominate the production of the HI only if $\zeta \gtrsim 4.5\times 10^{-16} \times (M_{\rm GMC}/10^6 \ M_{\odot})^{-1/2}$~s$^{-1}$, where $M_{\rm GMC}$ is the GMC mass, and including attenuation of the cosmic-ray fluxes. For most Galactic GMCs and conditions, FUV photodissociation dominates over cosmic-ray ionization for the production of the HI column densities. Furthermore, the cosmic-rays do not affect the HI-to-H$_2$ transition points.


INTRODUCTION
The compression of diffuse and warm interstellar atomic hydrogen (HI) gas into dense cold giant molecular (H 2 ) clouds (GMCs) is associated with radiative cooling, gravitational collapse, chemical complexity, and galaxy-star-and planet-formation across cosmic time (McKee & Ostriker 2007;Tacconi et al. 2020;Chevance et al. 2022;Sternberg et al. 2014, hereafter S14).Much of the cold (≲ 500 K) HI observed via 21 cm emissions and absorptions in the interstellar medium (ISM) of galaxies is produced in photodissociation regions (PDRs) in the atomic to molecular (HI-to-H 2 ) transition layers of the dense star-forming molecular clouds (Allen et al. 1986;Heiner et al. 2011;Walter et al. 2008;Bialy et al. 2017;Schruba et al. 2018;Saintonge & Catinella 2022).
What are the relative contributions of FUV photodissociation and cosmic-ray bombardment to the production of HI in typical molecular clouds in star-forming galaxies?In S14 and Bialy & Sternberg (2016, hereafter BS16) we presented numerical and analytic theory for the HI column densities produced by photodissociation in the HI-to-H 2 transition layers in optically thick and dusty PDRs, but with the exclusion of cosmicrays.In Bialy & Sternberg (2015) we investigated the HI/H 2 balance in low-metallicity cloud interiors dominated by cosmic-ray processes, but with no FUV.In Sternberg et al. (2021) we did consider combined FUV and cosmic-ray irradiation, but for dust-free systems in which cosmic-ray ionization, rather than dust catalysis, drives a gas-phase conversion of HI to H 2 , and in which the attenuation of the photodissocation rate is via pure H 2 absorption line self-shielding.Such dust-free PDRs may be relevant for young Universe conditions at the epoch of first star-formation.In this paper, we extend the analytic theory we presented in S14 and BS16 for dusty clouds, to also include cosmic-ray removal of the H 2 and the associated production of residual HI in the extended molecular cloud interiors.This in addition to direct photodissociation in the cloud surface PDRs.
In §2.1 and 2.2 we write down our basic HI/H 2 formation-destruction equation that includes a term for cosmic-ray removal of H 2 and the associated production of HI.We define the basic physical quantities and dimensionless parameters in the problem, α, β, σg , and G.We derive analytic expressions for the growth of the HI column density, from the outer PDR into the shielded cosmic-ray zone (CRZ), as a function of the gas density, far-UV field intensity, cosmic-ray ionization rate, and dust-to-gas ratio.In §2.3 we develop a formula for the critical cloud depths at which cosmic-rays dominate the the HI columns.In §2.4 we apply our formula to Galactic giant molecular clouds (GMCs) to assess whether cosmic rays can be significant contributors to HI columns in GMCs including their PDRs.In §3 we present numerical computations for the HI and H 2 abundance profiles and HI columns densities for a wide range of parameter combinations of the FUV intensity, cosmic-ray ionization rate, and gas density.We present results with and without the inclusion of a model for attenuation of the cosmic-ray fluxes.We also show how the profiles scale with the assumed dust-to-gas ratio.We discuss the effects of cosmic-ray ionization on the locations of the HI-to-H 2 transition points in the Appendix.We summarize in §4.

Formation-Destruction Equation
We consider an idealized one-dimensional semi-infinite cloud in slab geometry exposed on one side to beamed (normally incident) far-ultraviolet radiation, in combination with a flux of penetrating cosmic ray particles.In dusty systems the formation-destruction equation for the steady-state HI and H 2 fractions at any cloud depth is where x HI ≡ n HI /n is the atomic (HI) fraction, x H2 ≡ n H2 /n is the molecular (H 2 ) fraction, and n HI , n H2 , and n, are the atomic, molecular, and total hydrogen gas densities (cm −3 ).Particle conservation is where we assume that the abundances of hydrogen species other than HI or H 2 are negligibly small1 .The lefthand side of Eq. ( 1) is the H 2 formation rate (s −1 ), where is the grain-surface H 2 formation rate coefficient, T 2 ≡ T /(100 K) where T is the gas temperature (K), and σg is the dust-to-gas ratio normalized to the standard Galactic ISM dust-to-gas mass ratio of 1:100 for which σg =1 (Bohlin et al. 1978;Rémy-Ruyer et al. 2014).
The righthand side of Eq. ( 1), is the H 2 destruction rate by Lyman-Werner band photodissociation (LW: 912-1108 Å) in the PDR, and by cosmic-ray impact in the CRZ.In the first term, is the unattenuated free-space rate (Sternberg et al. 2014;Heays et al. 2017) for LW band photodissociation, where I UV is the far ultraviolet (6-13.6 eV) intensity relative to the Draine (1978) representation for the interstellar radiation field in the Solar neighborhood (I UV = 1).The factor of 1/2 accounts for the reduction of the photodissociation rate at the cloud surface due to the presence of the optically thick slab itself.The FUV and photodissociation rate are attenuated by a combination of dust absorption, and H 2 self-shielding as the LW absorption lines become optically thick.The exponential term in Eq. ( 1) accounts for the dust attenuation.The LW band dust optical depth where N = N HI + 2N H2 is the total (atomic plus molecular) hydrogen column density from the cloud surface, and is the dust absorption cross section per hydrogen nucleus.Here σg is the same dust-to-gas ratio appearing in Eq. ( 3).This parameter can also be viewed as the normalized dust absorption cross section.I.e., we are assuming that the H 2 formation rate coefficient (Eq.[3]) and the dust absorption cross section (Eq.[6]) scale identically with the overall dust abundance.
The H 2 self-shielding function where W d (N H2 ) (Hz) is the multi-line curve of growth for the H 2 dissociation bandwidth, and σ d = 2.36 × 10 −3 cm 2 Hz is the total H 2 dissociation cross section (see S14 for a detailed discussion of these quantities).We use the Draine & Bertoldi (1996) formula, as verified by S14, for the self-shielding function.At the cloud surface, N H2 = 0 and f = 1.For N H2 ≳ 10 14 cm −2 , the Doppler cores become optically thick and f becomes small.For N H2 ≳ 10 22 cm −2 the Lorentzian wings of the LW absorption lines overlap, andf → 0.
In the second term in Eq. (1), ϕζs(N ) is the local destruction rate of the H 2 by the cosmic rays.Here, is the unattenuated free-space rate of H 2 ionization by cosmic-ray impact This includes ionization by the primary cosmic-rays and the secondary energetic electrons.The parameter ϕ, of order unity, is the number of H 2 destruction events per cosmic-ray ionization.The H 2 destruction processes include ion-molecule chemical reactions driven by the initiating cosmic-ray ionizations, as well as direct cosmic-ray dissociation of the H 2 .In a steady-state, the hydrogen gas is primarily a mixture of HI and H 2 , and for a predominantly molecular medium ϕ ≈ 2 (Bialy & Sternberg 2015;Sternberg et al. 2021).
The factor C in the second term accounts for possibly different H 2 formation rates in the CRZs compared to the PDRs due to density and temperature gradients, as well as additional gas clumping in the CRZs.Differing formation rates imply For pressure equilibrium this then gives, For example, C ≈ √ 5 for pressure equilibrium between a FUV heated PDR with T ≈ 100 K, and a cosmic-ray heated CRZ with T ≈ 20 K. C can be increased further if there is any gas clumping.
The function s(N ) accounts for the possible attenuation of the cosmic-ray energy densities with cloud depth, and reduction of the associated H 2 ionization rates (Neufeld & Wolfire 2017;Padovani et al. 2018;Sternberg et al. 2021).However, the intrinsic energy spectra of the low-energy cosmic rays are uncertain, as are the transport mechanisms, e.g.free-streaming along magnetic field lines, or diffusive pitch-angle scattering off of pre-existing or self-generated MHD waves (Zweibel 2013;Padovani et al. 2020;Kempski & Quataert 2022).In our computations we either exclude cosmic-ray attenuation entirely (and set s = 1) or adopt a simple representative form for the cosmic-ray attenuation function.As in Sternberg et al. (2021), when including cosmic-ray attenuation we adopt the Padovani et al. (2018) broken power-law model Here N eff ≡ N/cosθ is the effective absorbing gas column density, where θ is the angle of the magnetic field along which the cosmic-rays propagate relative to the cloud normal, and N cr is the attenuation scale column.We use "model H" of Padovani et al. (2018) for which a = 0.385, and N cr = 10 19 cm −2 , and set cosθ = 1.This model is in agreement with observed declines of the cosmic-ray ionization rates with increasing cloud column densities (Caselli et al. 1998;Indriolo & McCall 2012;Neufeld & Wolfire 2017).See Fig. C1 in Padovani et al. (2022) for the full observational compilation.The power-law in Eq. ( 12) is valid for N eff between 10 19 and 10 24 cm −2 .For N eff < 10 19 cm −2 , s = 1, and divergence is avoided at small columns.Our basic question is: when (if ever) does internal cosmic-ray production of HI compete with photodissociation in the build-up of HI column densities in interstellar clouds?We are particularly interested in optically thick clouds consisting of fully developed outer photodissociation regions (PDRs) surrounding inner cosmic-ray dominated zones (CRZs).We are interested in the total HI columns, irrespective of the temperature-and depthdependent line widths of the associated 21 cm signatures.

ODE
For our 1D geometry, the atomic to molecular density ratio x HI /x H2 ≡ dN HI /dN H2 , and Eq. ( 1) can be written as the ordinary differential equation (ODE) In this equation the independent variable is N H2 (with N = N HI +2N H2 ) and the initial condition is N HI (0) = 0.The parameters and where n 2 ≡ n/(100 cm −3 ).Here and henceforth we assume ϕ = 2.The parameter α is the ratio of the unshielded free-space H 2 photodissociation rate to the molecular formation rate, and β is the ratio of the unattenuated cosmic-ray destruction rate to the molecular formation rate.
For characteristic interstellar conditions I UV ≈ 1, and α ≈ 1.9 × 10 4 for a cold gas density n 2 = 1.The parameter α remains large even for n 2 ≫ 1, especially near regions of active star-formation where the FUV field intensity I UV ≫ 1.At the cloud edge the H 2 is almost fully dissociated and the molecular fraction The molecular fraction grows as the FUV is attenuated with increasing cloud depth.The Galactic cosmic-ray ionization rate also varies depending on location, with values approaching 10 −15 s −1 in diffuse gas down to ∼ 10 −17 s −1 in dense clouds (e.g.Caselli et al. 1998;Indriolo & McCall 2012).See also the observational summary in Fig. C1 in Padovani et al. (2022).Some of this variation may be indicative of attenuation of the cosmic-ray fluxes as they traverse the clouds (Neufeld & Wolfire 2017).Here we adopt ζ −16 = 1 as a global characteristic value for the Galactic free-space cosmic-ray ionization rate.With ζ −16 = n 2 = T 2 = σg = 1, and with C = √ 5 for the CRZ, β = 3.0 × 10 −2 .The cosmic-ray ionization rates may be substantially larger in clouds near supernova remnants (Indriolo et al. 2010;Ceccarelli et al. 2011;Schuppan et al. 2012).
In most ISM environments α ≫ 1 and β ≲ 1.For constant Rn (independent of cloud depth) α is the maximal atomic to molecular density ratio at the fully photodissociated cloud edge, and β is the minimal atomic to molecular ratio in the optically thick cosmic-ray dominated interior (in the absence of CR attenuation).Complete atomic to molecular transitions are expected as the clouds become sufficiently optically thick.Given the solution to Eq. ( 13) for N HI (N H2 ), and with N ≡ N HI (N H2 ) + 2N H2 we obtain profiles for the HI column, N HI , and the derivatives, x HI and x H2 , as functions of N .In §3 we present such profiles computed numerically, but we first discuss several analytic solutions, as follows.

No CR attenuation
In the absence of any CR attenuation, with s ≡ 1 everywhere, and for any β, the HI fraction in the CRZ is, The atomic fraction x HI = 1/2 for β = 2.For a predominantly molecular CRZ, i.e. for β ≪ 1, the residual HI density is (18) independent of the total cloud density at any point (see also Solomon & Werner 1971;Li & Goldsmith 2003).For example, for ζ −16 = σg = 1, and a CRZ temperature T 2,CRZ = 0.2, the HI density in the CRZ is n HI,CRZ = 7.4 cm −3 .
For β ≲ 1, and with s ≡ 1, an excellent approximate analytic solution to Eq. ( 13) is The first term on the right is the HI column built up by photodissociation, and the second term is the HI due to cosmic-rays, both as functions of the molecular column N H2 .In this expression, where W g (N H2 ; σ g ) is the (universal) H 2 -dust limited curve of growth for the LW dissociation bandwidth (see S14 for a detailed discussion).For any σ g , W g is a preexisting function of N H2 , independent of the cloud parameters I UV or n.We use the analytic form for W g given by BS16 (their Eq. [27]).When all of the LW radiation is absorbed the integral converges to a constant, Here, W g,tot (σ g ) is the total dust-limited dissociation bandwidth (Hz), and G is then the (dimensionless) average H 2 self-shielding factor within an H 2 -dust absorption column.The righthand side of Eq. ( 21) is our BS16 fitting formula for G based on the multi-line (Meudon) PDR model computations we presented in S14.
At cloud depths beyond which all of the LW radiation is absorbed, Eq. ( 19) becomes where the basic dimensionless parameter The subscript t refers to optically thick.The first term in Eq. ( 22) is the total (asymptotic) HI column density produced by just photodissociation in optically thick PDRs, This is the formula for the total HI column density for beamed FUV fields derived by S14 in the absence of cosmic rays (i.e., for β = 0).The second term in Eq. ( 22) is the additional HI column produced by the cosmicrays, and it grows arbitrarily large with N unless the cosmic-ray ionization rate is sufficiently attenuated.In this term we have used the relation N H2 = N/(2 + β) for the CRZ in replacing N H2 with N in Eq. ( 19).
Differentiation2 shows that Eq. ( 19) is a good (but formally approximate) solution to Eq. ( 13) so long as βσ g N H2 is everywhere negligible compared to either σ g N HI or 2σ g N H2 .The latter two quantities are the HI-dust and H 2 -dust optical depths associated with the HI and H 2 columns respectively.Thus, if β = 0, i.e. with no cosmic-rays, Eqs. ( 19) and ( 22) are exact.But Eq. ( 19) remains accurate, so long as β < 1.We verify this in §3 by integrating Eq. ( 13) numerically and comparing to our analytic formulae.

With CR attenuation
When CR attenuation is included, the HI fraction in the CRZ decreases with cloud depth as CR attenuation reduces the HI that is built up in the CRZs.
Our analytic approximation for N HI,t can be generalized for arbitrary cosmic-ray attenuation functions, s(N ), by making the replacement for the second term in Eq. ( 22).In the second line we have evaluated the integral assuming the attenuation function given by Eq. ( 12) with cosθ = 1.

Critical Cloud Columns
The two terms on the righthand side of Eq. ( 22) are equal at the critical gas column, N = N crit , at which the cosmic rays start to dominate the growth of the HI column.Thus, for unattenuated cosmic-rays (s ≡ 1), Multiplying through by σ g gives the critical dust opacity, which depends on just the two (dimensionless) parameters αG and β.At a gas column N (or dust opacity τ dust ), cosmic-rays dominate the production of the HI if N > N crit (or if τ dust > τ dust,crit ), otherwise photodissocation dominates.In the weak-field limit, αG/2 ≪ 1 (and assuming β ≲ 1), we have3 In this limit most of the HI produced by photodissociation is built up past the HI-to-H 2 transition point in gas that is primarily molecular (as shown in Fig. 2 in § 3, see also Appendix).When cosmic-rays are added the CRZs and PDRs overlap, and the critical column is a measure of the relative HI production efficiency by photodissociation versus cosmic-ray ionization in the molecular gas.The critical column is therefore proportional to I UV /ζ −16 , independent of the cloud density n and/or H 2 formation rate.In the strong-field limit, αG ≫ 1 (and again for β ≲ 1), (31) For strong fields the photodissociated HI columns are built up in a (self-limited) fully atomic outer layer, and are only weakly (logarithmically) dependent on I UV .The cosmic-ray contributions to the HI occur in the inner fully optically thick regions where the atomic fractions are proportional to β.The critical gas column is therefore proportional to the ratio of the H 2 formation rate to ionization rate, or to the density to ionization rate for a given temperature, and the logarithmic factor of order unity.
In both the weak-and strong-field limits the critical columns are independent of the dust-to-gas ratio σg .
The intermediate case, αG/2 ≈ 1, is also important, because for a narrow range around this UV to gas density ratio (I UV /n 2 ≈ 3) a two-phased (WNM/CNM) thermal equilibrium is possible for fully atomic (HI) gas (Wolfire et al. 2003;Krumholz et al. 2008;Bialy & Sternberg 2019, S14).The range for two-phased equilibria is αG ∼ 1 to 4, weakly dependent on metallicity.Starforming gas and associated HI in the Milky Way and other galaxies may be self-regulated to be in a multiphased state (Ostriker et al. 2010).For such systems, the critical column is then (32) Here we are assuming that heating in the fully dissociated HI layers is via FUV photoelectric emission with negligible energy input by the cosmic-rays so that the thermal phase structure for the HI is primarily dependent on αG, i.e. on the ratio I UV /n (Wolfire et al. 2003;Bialy & Sternberg 2019).

Giant Molecular Clouds
Our expressions for the critical cloud columns are for one-sided illumination by a beamed (normally incident) FUV field, e.g. for a cloud irradiated by a nearby hot star.In the ambient medium two-sided irradiation by the background FUV field is more appropriate, and the critical columns are then doubled for a given cosmic-ray ionization rate.Furthermore, the irradiation may be isotropic rather than beamed.
For example, for a Galactic giant molecular cloud (GMC) embedded within ambient photodissociated HI containing a two-phased mixture of CNM and WNM the two-sided intermediate case αG = 2 applies.Galactic GMCs have characteristic hydrogen column densities N GMC ∼ 1.5 × 10 22 cm −2 (Solomon et al. 1987;McKee & Ostriker 2007;Lada & Dame 2020;Chevance et al. 2022) that are only weakly dependent on the cloud mass, M GMC , over a large range (∼ 10 to near 10 7 M ⊙ ) implying a mass-radius relation that scales approximately as M ∼ R 2 (Larson 1981).
For spherical systems irradiation by ambient isotropic FUV fields is the more natural configuration (e.g., Mc-Kee & Krumholz 2010).As discussed by S14 the HI column produced by photodissociation on one side of a plane-parallel slab exposed to isotropic 4 or a given H 2 photodissociation rate at the cloud surface, the incident radiation flux for isotropic fields is equal to half that for beamed fields, and αG is divided by 4 rather than 2 in Eq. ( 36).The factor ⟨µ⟩ = 0.8 is an angular average.See S14 for a detailed discussion.radiation is where in this expression ⟨µ⟩ = 0.8.(The subscript i refers to isotropic.)For a sphere, for which the PDR is a thin shell surrounding a CRZ core, the mean PDR HI column is (37) The mean CRZ HI column is where ⟨N ⟩ is the mean cloud column density.The critical mean column for a spherical cloud in an isotropic FUV field, and without CR attenuation, is then given by This expression is analogous to Eq. ( 28) with the extra factor of 2 for two sided beamed illumination of a slab.The critical β is hardly altered in switching from twosided beamed slabs to isotropically illuminated spheres.For example, for αG = 2, the critical β increases by just ∼ 10% for spheres.Eqs. ( 33) and ( 34) are therefore unaltered for spheres, but with N 22,GMC understood as the mean GMC column, as in Eq. ( 35).Setting n2,GMC T 2 C in Eq. ( 34), and with Eq. ( 35), we obtain the critical mass below which the GMCs are FUV dominated, and above which they are cosmic-ray dominated.With any gas clumping inside the GMC, or with the inclusion of CR attenuation, the critical masses will be larger still.

F
It follows from Eq. ( 40) that for a GMC temperature T 2,GMC = 0.2 the critical ionization rate scales as for a standard N 22,GMC = 1.5.More massive GMCs require lower ionization rates to be cosmic-ray dominated because their mean densities are lower.For example, for 6 × 10 6 M ⊙ near to the upper end of the Galactic GMC mass distribution (Williams & McKee 1997) ζ −16,crit ≈ 2 × 10 −17 s −1 .For a more typical GMC mass of 10 4 M ⊙ , the critical ionization rate is 5 × 10 −16 s −1 .
Relaxing the multiphase requirement, and deriving instead the critical column using Eq. ( 39) for weak (αG ≪ 1) isotropic FUV fields, GMCs are critical for or for independent of σg , and independent of the gas density or cloud mass.We stress again that Eqs. ( 42) and ( 43) hold for either slabs exposed to two-sided beamed fields or spheres illuminated by isotropic radiation.In Eq. ( 43) we have assumed that C = 1 since the PDRs and CRZs overlap in this limit (see § 3.1.1).Remarkably, in the weak-field limit, and for an ambient I UV ≈ 1, the critical ionization rate for typical GMCs is close to the characteristic Galactic ionization rate ζ −16 ≈ 1.Our analysis of the GMCs thus far does not include the effects of CR attenuation, which we do consider in § 3.3 below.CR attenuation increases the critical ionization rates further.

COMPUTATIONS
We now present numerical computations of the HI and H 2 density profiles and integrated HI column densities produced in gas slabs that are irradiated by combined fluxes of FUV photons and cosmic-rays.As is indicated by our analytic expressions Eqs. ( 19) and ( 22), the basic dimensionless parameter for the FUV driven HI-to-H 2 density profiles is αG (Eq.[23]) rather than α alone (see also S14).
For the cosmic-rays the basic parameter is β (Eq.[15]).In this paper we are focussing on the regime β ≲ 1 for which the gas is molecular in the absence of FUV, even without any CR attenuation.But the residual atomic component produced by the cosmic-ray bombardment contributes to the build up of the HI column densities.
We consider a wide range of conditions i.e., a range of αG, and β, for varying dust-to-gas abundance ratios σg , and we present results with and without the inclusion of cosmic-ray attenuation.We use Scipy ODEINT to integrate Eq. ( 13) and solve for x HI /x H2 as a function of gas column N ≡ N HI + 2N H2 , subject to x HI + 2x H2 = 1, for any αG, β and σg .When including CR attenuation we use the simple power-law form Eq. ( 12) for s(N ), assuming a normal magnetic field orientation cosθ = 1 in the definition of N eff .We compare our numerical integrations to our analytic approximations for N HI (N ) given by Eqs. ( 19) and ( 22).
We plot x HI and 2x H2 (blue and orange curves) as functions of the hydrogen column density N .The corresponding dust opacity, τ dust , is shown along the auxiliary (upper) x-axes.The black dashed curves are the depthdependent HI column densities found in our numerical integration of Eq. ( 13).The column density scale is shown along the righthand auxiliary y-axis.The overlying red dashed curves show the analytically computed HI columns.The agreement between the numerical and the analytically computed HI columns is excellent.
As expected, the hydrogen is primarily atomic at the cloud edges, and the molecular fractions are very small, with x H2 ≈ x H2 /x HI = 2/α = 6.2 × 10 −5 .Within the CRZs the gas is primarily molecular.In the absence of CR attenuation the HI fractions approach a cosmic-ray floor x HI = β/2 = 5 × 10 −2 (upper panel).
The red marker dots indicate the numerically computed HI-to-H 2 transition points, defined as the cloud depths where x HI = 2x H2 (or x HI = 1/2).For the models in Fig. 1 these occur at N tran = 1.4 × 10 20 cm −2 , or τ dust,tran = 0.26.The vertical dashed green lines indicate these positions as given by the analytic BS16 formula for the transition point (their Eq. [39]).We discuss this formula in the Appendix (Eq.[A1)]) and its continued range of applicability when cosmic-rays are included.
The vertical red dashed lines mark the columns, N 90 , where 90% of the incident FUV radiation is absorbed.This occurs at τ dust ∼ 1, and defines the inner edge of the PDR.The HI column produced by photodissociation is N HI,PDR = 2.1×10 20 cm −2 (see Eq. [24]).The vertical blue dashed lines mark the critical columns, N crit where the cosmic-rays start dominating the growth of the HI columns.Without CR attenuation this occurs at N crit = 4.4×10 21 cm −2 .The corresponding critical dust opacity is τ dust,crit = 8.4.
When CR attenuation is included (lower panel) the atomic fraction falls below the β = 0.1 cosmic-ray floor of 5 × 10 −2 without attenuation.The "knee" in the HI profile near N = 2 × 10 21 cm −2 is where the more slowly attenuating CR ionization processes take over from exponentially reduced photodissociation in producing the HI.The atomic fraction continues to decline at greater cloud depths as βs(N )/2, and becomes very small.The much reduced HI abundance in the CRZ when CR attenuation is included moderates the growth of the HI column density (see Fig. 1) and the critical column is now N crit = 7.1 × 10 22 cm −2 , or τ dust,crit = 135.3.When CR attenuation is included the CRZ must be 135 times larger than the PDR for cosmic-rays to contribute signifcantly to the production of the HI.In Fig. 2 we present an αG versus β model grid for the HI-to-H 2 density profiles, and integrated HI column densities, assuming σg = 1, and no attenuation of the cosmic-ray ionization rates (s = 1).From top to bottom, αG ranges from 0.01 (weak-field limit) to 10 (strong field limit).From left to right β ranges from 0 to 1, i.e., weak to moderate5 cosmic-ray irradiation.For σg = 1 these ranges correspond to I UV /n 2 from 1.7 × 10 −2 to 17.0, and ζ −16 /n 2 from 0 to 33.3 (for C = √ 5, T 2 = 1, and ϕ = 2).
As in Fig. 1, in each panel we show the HI and H 2 fractions, x HI (blue curves) and 2x H2 (orange curves), as computed by integrating Eq. ( 13) numerically.The black and red dashed curves are the HI column densities found in our numerical integrations and using our analytic formulae respectively.The agreement between the two curves is excellent across the entire parameter space.
Again, the hydrogen is atomic at the cloud edges, and the molecular fractions are very small, with x H2 ≈ x H2 /x HI = 2/α from 6.2 × 10 −3 to 6.2 × 10 −6 .Within the CRZs the gas is (by assumption) primarily molecular, and the HI fractions approach x HI = β/2, i.e. range from 5 × 10 −4 to 0.5, from small to moderate cosmic-ray ionization rates.The red dots are the HI-to-H 2 transition points.For the range of αG in Fig. 2 these occur at gas columns, N tran , equal to 1.9 × 10 17 , 3.9 × 10 18 , 1.4 × 10 20 , and 9.1 × 10 20 cm −2 , corresponding to dust optical depths, τ dust , equal to 3.6 × 10 −4 , 7.4 × 10 −3 , 0.26, and 1.7.The vertical green dashed lines show these positions using the BS16 formula, Eq. (A1).In the weak-field limit (small αG) an HI-to-H 2 transition is induced by H 2 self-shielding at small cloud depths where τ dust ≪ 1, and dust attenuation is irrelevant for the transition point.Most of the photodissociated HI column density is built up inside the predominantly molecular zone, up to τ dust ≈ 1 where the FUV is finally fully absorbed.In the strong field limit (large αG) the fully atomic layer becomes sufficiently large that the dust associated with this layer (the "HI-dust") dominates the absorption of the FUV.The transition to H 2 is then very sharp, and most of the HI column is produced in the outer fully dissociated layer.Because we are assuming β ≤ 1 the cosmic-rays do not inhibit transitions to molecular gas as the clouds become optically thick to the FUV, and the transition points are unaffected (see Appendix).
The vertical red dashed lines shown for the β = 0 cases (leftmost column in Fig. 2) show the FUV "absorption columns", N 90 , where 90% of the photodissociated HI columns are built up.The 90% absorption depths occur at τ dust ≈ 1, independent of αG, and are unaffected by the presence of cosmic-rays.We do not display the N 90 lines for the β ̸ = 0 panels.Instead, for β ̸ = 0 the vertical blue dashed lines indicate the critical gas columns, N crit , and dust opacities, τ dust,crit , where the cosmic-ray and FUV contributions to the integrated HI columns are equal.
In Fig. 3, we plot curves as given analytically by Eq. ( 29) for the critical dust opacities, τ dust,crit , as functions of αG and β.The blue squares are the critical opacities as found numerically in Fig. 2. They lie very close to the analytic curves.The auxiliary y-axes in Fig. 3 show the corresponding critical gas column densi- HI-to-H2 density profiles as functions of the gas column N (lower x-axes) and the dust optical depth τ dust (upper x-axes) for αG from 0.01 (weak field) to 10 (strong field), and for cosmic-ray parameters β from 0 to 1, and with no cosmic-ray attenuation.The gas-to-dust ratio σg = 1.The curves are for the HI fractions xHI (blue), twice the H2 fraction 2xH 2 (orange), and the HI column density NHI, integrated numerically (dashed black), and using our analytic formula Eq. ( 19) (dashed red).The red dots mark the HI-to-H2 where xHI = 2xH 2 .The analytic approximation for the transition points (Eq.[A1]) are indicated by the vertical dashed green lines.For β = 0 (left column) the vertical dashed red lines are for N90 where 90% of the photodissociated HI columns are built up.For β ̸ = 0 the vertical dashed blue lines mark the critical cloud depths where the cosmic-ray contributions to the HI columns are equal to the photodissociated HI columns.
ties assuming σg = 1 (Eq.[28]).The left panel displays τ dust,crit as a function of β (or ζ −16 /n 2 C for σg = 1) for several values of αG from 0.01 to 100.The right panel displays curves for τ dust,crit as a function of αG (or I UV /n 2 for σg = 1) for several values of β from 0.001 to 1.The curves illustrate the limiting behaviors given by Eqs. ( 30) and (31).For a given αG the critical dust opacities and gas columns always vary inversely with β.For a given β, they vary linearly with αG in the weakfield limit (αG ≪ 1) and logarithmically with αG in the strong-field limit (αG ≫ 1).
The horizontal dashed blue line in Fig. 3 is the τ dust = 1 boundary between the PDR and the CRZ.The curves again show that in the weak-field limit, αG ≪ 1 cosmicray production of the HI can become competitive with photodissociation already within the PDRs (i.e.within τ dust ≲ 1).Conversely, in the strong-field limit, αG ≫ 1, the critical opacities become large with τ dust,crit > 1, even if β approaches 1.In this limit the cosmic-ray production of the HI occurs mainly in the optically thick cloud interiors.

HI and H2 Profiles
In Fig. 4 we display the same αG versus β grid for the HI and H 2 profiles as in Fig. 2 but now with the inclusion of cosmic-ray attenuation.We again assume σg = 1.In all panels, we assume the broken power-law CR attenuation function s(N ) as given by Eq. ( 12), with cosθ = 1 for the magnetic field orientation.The effect of the cosmic-ray attenuation is most clearly seen for β = 1 in the righthand column of Fig. 4. Without attenuation the HI fraction x HI = 1/3 at large depths for β = 1 (see Fig. 2), and the integrated HI columns therefore rise sharply with increasing cloud depth.With attenuation the local HI fractions decrease and the resulting integrated HI columns are reduced.
The black dashed curves in Fig. 4 show the HI columns found by numerically integrating Eq. ( 13) (again using Figure 3. Critical dust opacities, τ dust,crit , at which the cosmic-ray contributions to the HI column densities are equal to the photodissociated HI columns.Right panel: τ dust,crit as functions of αG with individual curves as given by our analytic expression Eq.( 29), for β from 0.001 to 1.The blue squares are the numerically computed critical depths shown in Fig. 2. The auxiliary axes for IUV/n2 and Ncrit are for a dust-to-gas ratio σg = 1.The horizontal blue dotted line is for τ dust = 1 below which photodissociation and cosmic-ray production of the HI overlap (see text).The horizontal red dotted line corresponds to the typical half-column density of Galactic GMCs.The vertical green line marks the intermediate αG ≈ 2 regime for which multiphased HI is possible, and within the grey strip for which multi-phased HI is possible for FUV heated gas.Left panel: τ dust,crit as functions of β with individual curves as given by our analytic expression Eq.( 29), for αG from 0.01 to 100.The blue squares are the numerically computed critical depths shown in Fig. 2. The auxiliary axes for ζ−16/n2C and Ncrit are for a dust-to-gas ratio σg = 1.The horizontal blue dotted line is for τ dust = 1 below which photodissociation and cosmic-ray ionization overlap.The horizontal red dotted line corresponds to the typical half-column density of Galactic GMCs.
Scipy ODEINT), but now including the attenuation function s(N ).The red dashed curves show the HI columns computed using our analytic approximation Eq. ( 19) using Eq. ( 27) for the CR term.The agreement between the numerical solution and the analytic representation is excellent.
As in Fig. 2 the vertical blue lines in Fig. 4 mark the critical cloud depths at which the FUV and CR contributions to the HI columns are equal.Due to the reductions in the HI fractions in the CRZs the critical depths are increased compared to the no CR attenuation case.The effect is especially significant in the strong FUV field limit αG > 1 for which the FUV contributions to the HI columns become large.(In some of the panels in Fig. 2 the blue markers do not appear because the critical depths are off scale high).
The red dots in Fig. 4 show the HI-to-H 2 transition points where x HI = 2x H2 .The vertical green lines mark the transition points as estimated using the BS16 formula Eq. (A1).The positions of the transition points are fully controlled by the FUV radiation absorption, and are not affected by the presence of cosmic rays or the inclusion of CR attenuation.

Critical Dust Opacities and Gas Columns:
With CR Attenuation In Fig. 5, we plot curves for the critical dust opacities, τ dust,crit , as functions of αG and β, but now with the inclusion of CR attenuation as for the profiles shown in Fig. 4. To generate these curves we modify Eq. ( 29) for τ dust,crit by making the replacement given by Eq. ( 27) in Eq. ( 22).The left panel displays curves for τ dust,crit as functions of β (or ζ −16 /n 2 C for σg = 1) for several values of αG from 0.01 to 100.The right panel shows τ dust,crit versus αG (or I UV /n 2 for σg = 1) for several values of β from 0.001 to 1.The auxiliary y-axes in Fig. 3 show the corresponding critical gas column densities assuming σg = 1.The blue squares are the results of the numerical integrations found in Fig. 4, and they lie very close to the analytic curves.The primary affect of CR attenuation is to steepen the critical curves, since attenuation dampens the growth of the HI columns preferentially at low β and large αG.
For example, for αG = 0.01 and β = 0.001, τ dust,crit is increased to from 10 to 200 when CR attenuation is included, with N crit increasing to 1.1 × 10 23 cm −2 .As β is increased for αG = 0.01, the critical points move  inward, the PDRs and CRZs overlap as seen in the top row of Fig. 4, and the attenuation effects are reduced due to the rapid build up of the CR contributions.As another example, for αG = 1 and β = 0.1, τ dust,crit increases from 8.5 to 135. with N crit increasing to 7.1 × 10 22 cm −2 .

GMCs and Multiphased HI
The horizontal red dashed lines in Figs. 3 and 5 mark the half-column, N GMC /2 = 7.5 × 10 21 cm −2 , for typical Galactic GMCs, as discussed in § 2.4.For any αG and β for which N crit < N GMC /2, the GMC is "supercritical" and the CRZ dominates the total HI column density.For N crit > N GMC /2 the GMCs are "subcritical" and the PDR dominates the HI.GMCs are just critical for αG and β at the intersections of the critical curves with the N GMC /2 line.
As discussed in § 2.4, without CR attenuation and in the weak-field limit GMCs are critical for β ≈ 7.0 × 10 −2 αG, or for ζ −16 ≈ 0.6I UV (Eqs.[42] and [43]).This relation is seen in Fig. 3 moving along the red line for small αG.As shown in Fig. 5, with CR attenuation the critical β and ζ −16 are much larger.For example, for αG = 0.1, and without CR attenuation β = 7.0 × 10 −3 for critical GMCs, and this increases to 6.0 × 10 −2 when CR attenuation is included.Or, for I UV = 1, the critical ionization rate ζ −16 increases from 0.6 to 5.1 for models with and without CR attenuation respectively.
The vertical green dashed lines in the righthand panels of Fig. 3 and 5 mark the intermediate αG = 2 case (nominally I UV /n 2 ≈ 3), for which multiphased (WNM/CNM) HI is possible in the PDRs, as indicated by the grey strip.Without CR attenuation, and at αG = 2, the red GMC line in Fig. 3 intersects the critical curve for β = 0.1.This is as given by Eq. (33).Fig. 5 shows that β = 0.9 when CR attenuation is included for critical GMCs with αG = 2.For example, for the nominal I UV /n 2 ≈ 3, the critical ratio ζ −16 /(n 2 C) increases from 1.5 to 13.5.The critical free-space ionization rate then scales with GMC mass as when CR attenuation is included.
As discussed in §2.4, the critical columns and ionization rates for spheres illuminated isotropically are essentially identical to the critical values for two-sided slabs, where the gas column N for slabs is replaced by the mean column ⟨N ⟩ for spheres.This is because 2N HI,PDR,i ≈ N HI,PDR (see Eqs. [24] and [36]).
In Fig. 6 we further consider the αG = 2 case.We show the HI fractions x HI (blue dashed curves), and integrated HI column densities N HI (black curves), for several values of β.In the upper panel CR attenuation is excluded, and in the lower panel CR attenuation is included.For both the dust-to-gas ratio is σg = 1.For αG = 2 the HI column produced by photodissociation is N HI,FUV = 3.65 × 10 20 cm −2 .The horizontal green lines are at twice this value (see the righthand column density scale) and the intersections with the N HI curves are at the critical cloud depths for each β.The vertical red lines mark the typical half-column of 7.5×10 20 cm −2 for the Galactic GMCs.Without CR attenuation we again see that GMCs are critical for β = 0.1.With CR attenuation the critical value is much larger with β = 0.9.This corresponds to a very large free-space ionization rate to density ratio ζ −16 /(n 2 CT 1/2 2 ) = 13.4.

Dust-to-Gas Ratio
How do the HI-to-H 2 profiles depend on the assumed dust-to-gas ratio, as parameterized by our σg ?The dustto-gas ratio (as controlled by the overall metallicity) enters in two ways.First via the H 2 dust-grain formation rate coefficient R (eq.[3]), and second via the FUV dust absorption cross section σ g (Eq.[6]).The formation rate coefficient, R, appears in the denominators of our dimensionless parameters α and β (Eqs.[14] and [15]) in our ODE Eq. ( 13).The dust absorption cross section, σ g , appears in the definition of the dust optical depth, τ dust , in Eq. ( 13).But importantly, the fundamental parameter αG is only weakly dependent 6 on σg due to the cancellation when taking the ratio σ g /R (see Eq. [23]).
BS16 studied the β = 0 case (i.e.no cosmic rays) and found that when expressed in terms of τ dust (rather than the gas column N ) then to a very good approximation, especially for σg in the range 0.1 to 10, the HI-to-H 2 transition points depend on just αG independent of σg .This is the essence of our Eq.(A1).This invariance is somewhat surprising especially in the weak-field limit where the FUV attenuation is governed purely by H 2 self-shielding, τ dust,tran ≪ 1, and dust-shielding plays no role.Furthermore, at cloud depths beyond the transition points, i.e. within the molecular zones, the HI and H 2 density profiles as functions of τ dust , also depend on just αG, and are invariant with σg to a very good approximation.Within the fully atomic outer layers, and up to the invariant transition points, the H 2 profiles do depend on σg , with H 2 fractions at the optically thin cloud edges that vary inversely with the dust abundance and the associated H 2 formation rate coefficient.Fig. 7 illustrates the behavior when cosmic rays are included.In this example we set αG = 1 and β = 0.1 and present results for σg = 0.1, 1, and 10.In the upper panel we exclude CR attenuation, and in the lower panel CR is included according to Eq. ( 12).Once again, the blue and orange curves are the atomic and molecular fractions x HI and 2x H2 , and the black curves are the integrated HI column densities.The red markers indicate the transition points, as do the vertical green dashed lines according to Eq. (A1).
The cosmic-ray ionization rate is unaffected by the presence of dust, and in the absence of CR attenuation the HI profiles within the optically thick CRZs do not depend on σg and the atomic fraction reaches the cosmic-ray floor x HI = 4.8 × 10 −2 (see upper panel).However, with CR attenuation, and when expressed as functions of the dust optical depth, the HI fractions in-crease with σg for a given τ dust (see lower panel).This is simply because the CR attenuation depends on the gas column N = τ dust /σ g .For our assumed CR attenuation power-law (Eq.[12]) with a = 0.385, x HI varies as σ0.385 g at a given dust optical depth.In all cases the atomic fractions decrease with cloud depth ([Eq. 26]).
We have verified by explicit computation that this ovearall behavior is maintained for the entire range of αG and β in our model grids for σg from 0.1 to 10.

SUMMARY
In this paper we extend the analytic treatment presented by Sternberg et al. (2014) and Bialy & Sternberg (2016) (S14 and BS16) for the production of atomic hydrogen (HI) via FUV photodissociation at the boundaries of interstellar molecular (H 2 ) clouds, to also include the effects of penetrating (low-energy) cosmic-rays for the growth of the total HI column densities.
We focus on idealized one-dimensional gas slabs, consisting of outer photodissociation regions (PDRs) and inner cosmic-ray zones (CRZs).We compute the depth dependent steady-state abundances of the HI and H 2 , in a balance between grain-surface formation of the H 2 and destruction via FUV photodissociation and cosmicray ionization.The FUV photodissociation rates are reduced by (standard) H 2 self-shielding, and dust absorption.For the cosmic-rays we assume either constant overall ionization rates, or models that include depthdependent attenuation of the cosmic-ray fluxes.
The physical parameters in the problem are (a) the free-space intensity, I UV , of the FUV radiation and the associated H 2 photodissociation rate; (b) the free-space cosmic ray H 2 ionization rate, ζ; (c) the density, n, of hydrogen nuclei, in atoms and molecules; (d) the H 2 formation rate coefficient R; (e) the FUV dust absorption cross section σ g ; (f) the gas temperature T ; and (g) a density enhancement factor, C, for the cool CRZs relative to the warmer PDRs.An additional (chemical) parameter is the number, ϕ, of H 2 dissociations per cosmic-ray ionization event.
The governing HI/H 2 formation-destruction equation that we solve is Eq. ( 1), or in differential form Eq. ( 13).The solutions for the HI and H 2 density profiles and the integrated HI columns, depend primarily on the ratios I UV /Rn and ζ/Rn, as encapsulated in our dimensionless parameters αG, and β (Eqs.[14], [15] and [23]).A third dimensionless parameter is the dust-to-gas ratio σg .It sets the magnitude of both the dust absorption cross section, and the molecular formation rate coefficient.
We solve Eq. ( 13) numerically, and we also develop simple analytic formulae for the growth of the HI column density in terms of αG and β (Eqs.[19] and [22]).Our analytic formulae provide an excellent match to the numerical integrations.Our focus is on conditions (β ≤ 1) for which the gas is primarily molecular in the optically thick cloud interiors.As we discuss in the Appendix, for these conditions cosmic-rays do not affect the locations of the HI-to-H 2 transition points.We consider both weak fields (αG ≪ 1) and strong fields (αG ≫ 1), and compute the critical cloud columns, N crit , at which cosmic-rays dominate the production of the total HI columns.We write down analytic expressions for the critical columns.We also examine how the HI and H 2 profiles scale with the assumed dust-to-gas ratio.
As an example, we apply our theory to Galactic giant molecular clouds (GMCs), with typical hydrogen gas column densities ∼ 1.5 × 10 22 cm −2 (independent of mass).For GMCs we consider both plane-parallel slabs exposed to beamed FUV fields, and spherical clouds illuminated by isotropic radiation.For weak FUV fields, for which I UV /n ≪ 3.4 × 10 −2 cm 3 , and with I UV = 1, the CRZ dominates the production of the HI if the free-space ζ > 5.1 × 10 −16 s −1 .This estimate for the critical ionization rate includes cosmicray attenuation within the GMCs.For multiphased warm/cold HI within the PDRs, for which I UV /n ≈ 3.4 × 10 −2 cm 3 , the CRZ dominates the HI if ζ ≳ 4.5 × 10 −16 × (M GMC /10 6 M ⊙ ) −1/2 s −1 , where M GMC is the GMC mass.The very large critical ionization rates suggest that FUV photodissociation dominates the production of the HI in most Galactic GMCs.points.In particular, for characteristic Galactic conditions, with ζ −16 /I UV ≈ 1 and σg = 1, this never occurs.For example, for αG = 2, Eq.(A4) gives β max = 7, so that for I UV = 1, the free-space cosmic-ray ionization rate must be larger than 3 × 10 −15 s −1 for the transition point to be affected.Even without CR attenuation the ionization rate would have to be very large, and greater than 9 × 10 −16 s −1 .

Figure 1 .
Figure 1.HI and H2 density profiles, for αG = 1 = 0.59IUV/n2, with T2 = 1 and σg = 1 (see Eq.[23]) and β = 0.1 = 3.0 × 10 −2 ζ−16/n2, with C = √ 5 (see Eq. [15]).The upper panel are results without CR attenuation.The lower panel includes CR attenuation.The profiles are functions of the gas column N (lower x-axes ) and the dust optical depth τ dust (upper x-axes).The curves are for the HI fractions xHI (blue), twice the H2 fractions 2xH 2 (orange), and the HI column densities NHI integrated numerically (dashed black), and using our analytic formula Eq. (19) (dashed red).The left y-axes show the fractions, and the right y-axes are the column density scales.The red dots mark the HI-to-H2 transition points where xHI = 2xH 2 .The analytic approximations for the transition points (Eq.[A1]) are indicated by the vertical dashed green lines.The vertical dashed red line is N90 where 90% of the photodissociated HI columns are built up.The vertical dashed blue lines mark the critical cloud columns Ncrit where the cosmic-ray contributions to the HI columns start to dominate.

3. 1 .
Model Grid: No CR Attenuation, and σg = 1 3.1.1.HI and H2 Profiles Figure 2.HI-to-H2 density profiles as functions of the gas column N (lower x-axes) and the dust optical depth τ dust (upper x-axes) for αG from 0.01 (weak field) to 10 (strong field), and for cosmic-ray parameters β from 0 to 1, and with no cosmic-ray attenuation.The gas-to-dust ratio σg = 1.The curves are for the HI fractions xHI (blue), twice the H2 fraction 2xH 2 (orange), and the HI column density NHI, integrated numerically (dashed black), and using our analytic formula Eq. (19) (dashed red).The red dots mark the HI-to-H2 where xHI = 2xH 2 .The analytic approximation for the transition points (Eq.[A1]) are indicated by the vertical dashed green lines.For β = 0 (left column) the vertical dashed red lines are for N90 where 90% of the photodissociated HI columns are built up.For β ̸ = 0 the vertical dashed blue lines mark the critical cloud depths where the cosmic-ray contributions to the HI columns are equal to the photodissociated HI columns.

Figure 4 .
Figure 4. Model grid with curves and markers as in Fig. 2, but with the inclusion of CR attenuation.

Figure 5 .
Figure5.As in Fig.3for the critical dust opacities, but with the inclusion of CR attenuation.

Figure 6 .
Figure6.HI fractions xHI (dashed blue curves) and integrated HI columns NHI (black curves) versus cloud depth, as parameterized by the gas column, N , or the dust optical depth, τ dust , for the intermediate, beamed field, αG = 2 case for multiphased HI, for a dust-to-gas ratio σg = 1.The curves are labelled by the assumed values of β.In the upper panel CR attenuation is not included, and in the lower panel CR attenuation is included.The horizontal green line is for an HI column equal to twice the photodissociated column for αG = 2.The vertical red line marks the half gas column for typical Galactic GMCs.

Figure 7 .
Figure 7. HI and H2 fractions, xHI and 2xH 2 (blue and orange curves) and HI column densities, NHI (black curves) assuming αG = 1, and β = 0.1, for σg =0.1, 1, and 10, without and with CR attenuation (upper and lower panels).The curves are plotted as functions of the dust optical depth τ dust , or equivalently σg × N , where N is the gas column density.