A Simple and Accurate Network for Hydrogen and Carbon Chemistry in the Interstellar Medium

2.1.1. General Framework

The chemical network in this paper is listed in Tables 1 and 2, and described in detail in Appendix A. Our chemical network is based on the NL99 network in Nelson & Langer (1999) and Glover & Clark (2012a), with significant modification and extension. Eighteen species are considered: H, ${{\rm{H}}}_{2}$, ${{\rm{H}}}^{+}$, ${{\rm{H}}}_{2}^{+}$, ${{\rm{H}}}_{3}^{+}$, He, He⁺, O, ${{\rm{O}}}^{+}$, C, ${{\rm{C}}}^{+}$, CO, HCO⁺, Si, ${\mathrm{Si}}^{+}$, e, ${\mathrm{CH}}_{x}$, and ${\mathrm{OH}}_{x}.$ Following NL99, ${\mathrm{CH}}_{x}$ (including CH, CH₂, CH⁺, ${\mathrm{CH}}_{2}^{+}$, ${\mathrm{CH}}_{3}^{+}$) and ${\mathrm{OH}}_{x}$ (including OH, H₂O, OH⁺, ${{\rm{H}}}_{2}{{\rm{O}}}^{+}\,{{\rm{H}}}_{3}{{\rm{O}}}^{+}$) are pseudo-species, and their treatment is described in Appendices A.2 and A.3. Among the 18 species in our network, 12 of them are independently calculated, and 6 of them, H, He, C, O, Si and e, are derived from conservation of total hydrogen, helium, carbon, oxygen, and silicon nuclei, and the conservation of total charge.³ The gas-phase abundance of total helium nuclei is ${x}_{\mathrm{He},\mathrm{tot}}\,={n}_{\mathrm{He},\mathrm{tot}}/n=0.1$. The gas-phase abundances of total carbon, oxygen, and silicon nuclei are assumed to be proportional to the gas metallicity relative to the solar neighborhood ${Z}_{{\rm{g}}}$, and we adopt the values ${x}_{{\rm{C}},\mathrm{tot}}=1.6\times {10}^{-4}{Z}_{{\rm{g}}}$ (Sofia et al. 2004), ${x}_{{\rm{O}},\mathrm{tot}}=3.2\times {10}^{-4}{Z}_{{\rm{g}}}$ (Savage & Sembach 1996), and ${x}_{\mathrm{Si},\mathrm{tot}}=1.7\times {10}^{-6}{Z}_{{\rm{g}}}$ (Cardelli et al. 1994). We do not explicitly follow the evolution of dust grains, but instead assume that the grain-assisted reaction rates scale with the dust abundance relative to the solar neighborhood Z_d. We assume that the gas metallicity Z_g and dust abundance Z_d vary simultaneously, and use a single parameter for metallicity relative to the solar neighborhood,

$$\begin{eqnarray}&&Z={Z}_{{\rm{g}}}={Z}_{{\rm{d}}}.\end{eqnarray} \tag{ 1 }$$

We have updated all reaction rates in the original NL99 network, according to the most recent values in the UMIST (McElroy et al. 2013) and KIDA (Wakelam et al. 2010) catalogs, and other references in Tables 1 and 2. Notably, we have a higher rate of ${{\rm{C}}}^{+}+\mathrm{OH}\to {\mathrm{CO}}^{+}+{\rm{H}}$, and a higher rate of ${{\rm{C}}}^{+}+{\rm{e}}\to {\rm{C}}$ by including both radiative and dielectronic recombination. These higher rates lower the electron abundance and aid the formation of ${\mathrm{OH}}_{x}$ and HCO⁺, resulting in more efficient CO formation.

We have also modified and extended the NL99 network, as described in detail in Appendix A.1. One important extension is the addition of grain-assisted recombination of ${{\rm{C}}}^{+}$ and He⁺. Although Glover & Clark (2012a) concluded these grain reactions do not affect CO formation at the low cosmic-ray ionization rate ${\xi }_{{\rm{H}}}={10}^{-17}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$, we found that with the updated cosmic-ray ionization rate ${\xi }_{{\rm{H}}}=2\times {10}^{-16}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$ (Indriolo et al. 2007, 2015), these reactions are critical for CO formation at moderate densities $n\lesssim 1000\,{\mathrm{cm}}^{-3}$. This is because electrons formed by cosmic-ray ionization inhibit CO formation, while He⁺ ions formed by cosmic rays destroy CO (Bisbas et al. 2015). Grain-assisted recombination reduces the abundance of e and He⁺. With these extensions, we have 50 reactions in total, including collisional reactions, grain-assisted recombination and ${{\rm{H}}}_{2}$ formation, photodissociation by FUV radiation, and cosmic-ray ionization.

2.1.2. Photochemistry

The photodissociation reactions induced by FUV depend on the radiation field strength. We assume that the incident radiation field scales with the standard interstellar radiation field determined by Draine (1978), using the parameter χ as the field strength relative to ${J}_{\mathrm{FUV}}=2.7\times {10}^{-3}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ (G₀ = 1.7 in Habing 1968 units). The photochemistry reaction rates appropriate for this radiation field are listed in Table 2.

In optically thick regions, the radiation is attenuated by dust and by molecular line shielding. In a plane-parallel slab geometry with beamed incident radiation field from one direction, the photodissociation rates in optically thick regions ${R}_{\mathrm{thick}}$ can be related to the rates in optically thin regions ${R}_{\mathrm{thin}}$ by

$$\begin{eqnarray}&&{R}_{\mathrm{thick}}=\chi {R}_{\mathrm{thin}}{f}_{\mathrm{shield}}=\chi {R}_{\mathrm{thin}}{f}_{\mathrm{dust}}{f}_{{\rm{s}}},\end{eqnarray} \tag{ 2 }$$

where ${f}_{\mathrm{dust}}=\exp (-\gamma {A}_{V})$ is the dust-shielding factor and ${f}_{{\rm{s}}}$ is the self-shielding factor of C, CO, and ${{\rm{H}}}_{2}$ (see text below). The values of the parameter γ appropriate for the different reactions are listed in Table 2. The visual extinction A_V is calculated as

$$\begin{eqnarray}&&{A}_{V}=\displaystyle \frac{{{NZ}}_{{\rm{d}}}}{1.87\times {10}^{21}\,{\mathrm{cm}}^{-2}},\end{eqnarray} \tag{ 3 }$$

where $N={N}_{{\rm{H}}}+2{N}_{{{\rm{H}}}_{2}}$ is the total hydrogen column density along the line of sight. This corresponds to ${R}_{V}\,\equiv {A}_{V}/E(B-V)=3.1$ and $N/E(B-V)=5.8\times {10}^{21}\,{\mathrm{cm}}^{-2}$, appropriate for the diffuse ISM (Bohlin et al. 1978). In clouds with slab geometry and isotropic radiation field impinging from one side (Section 4), ${f}_{\mathrm{shield}}$ is an average value calculated from different incident angles. For more general, non-slab geometry, photo rates at a given location would be calculated by angle averages of Equation (2), as both χ and ${f}_{\mathrm{shield}}$ vary with direction.

For photoelectric heating on dust (see Section 2.2), the heating rate depends on the radiation field strength between $6\mbox{--}13.6\,\mathrm{eV}$, which affects the charge states of grains. For radiation in the FUV important for the photoelectric effect, we use a constant dust cross-section ${\sigma }_{{\rm{d}},\mathrm{PE}}={10}^{-21}{Z}_{{\rm{d}}}\,{\mathrm{cm}}^{2}\,{{\rm{H}}}^{-1}=1.87{A}_{V}/N$ to calculate the attenuation by dust:

$$\begin{eqnarray}&&{\chi }_{\mathrm{PE}}(N)=\chi {f}_{\mathrm{dust},\mathrm{PE}}=\chi \exp (-{\sigma }_{{\rm{d}},\mathrm{PE}}N);\end{eqnarray} \tag{ 4 }$$

our adopted cross-section is close to the value $1.8{A}_{V}/N$ used in previous PDR models.

In addition to dust shielding, we also include the ${{\rm{H}}}_{2}$ self-shielding, CO self-shielding, C self-shielding, and shielding by ${{\rm{H}}}_{2}$ of CO and C. The photodissociation rate of ${{\rm{H}}}_{2}$ can be written as

$$\begin{eqnarray}&&{R}_{\mathrm{thick},{{\rm{H}}}_{2}}=\chi {R}_{\mathrm{thin},{{\rm{H}}}_{2}}{f}_{\mathrm{dust}}{f}_{{\rm{s}},{{\rm{H}}}_{2}}({N}_{{{\rm{H}}}_{2}}).\end{eqnarray} \tag{ 5 }$$

We use the results in Draine & Bertoldi (1996) for H₂ self-shielding:

$$\begin{eqnarray}&&{f}_{s,{{\rm{H}}}_{2}}({N}_{{{\rm{H}}}_{2}})=\displaystyle \frac{0.965}{{(1+x/{b}_{5})}^{2}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&+\displaystyle \frac{0.035}{{(1+x)}^{0.5}}\exp [-8.5\times {10}^{-4}{(1+x)}^{0.5}],\end{eqnarray} \tag{ 7 }$$

where $x\equiv {N}_{{{\rm{H}}}_{2}}/(5\times {10}^{14}\,{\mathrm{cm}}^{-2})$ and ${b}_{5}=b/(\mathrm{km}\,{{\rm{s}}}^{-1})$ for b the velocity dispersion. In regions where $2x({{\rm{H}}}_{2})=2n({{\rm{H}}}_{2})/n\gtrsim 0.1$, the self-shielding of H₂ is in the wings of the line profile, and the H₂ fraction is not sensitive to the choice of b₅. Here we use a constant value ${b}_{5}=3.$

Similarly, the photodissociation rate of CO can be written as

$$\begin{eqnarray}&&{R}_{\mathrm{thick},\mathrm{CO}}=\chi {R}_{\mathrm{thin},\mathrm{CO}}{f}_{\mathrm{dust}}{f}_{{\rm{s}},\mathrm{CO}}({N}_{\mathrm{CO}},{N}_{{{\rm{H}}}_{2}}).\end{eqnarray} \tag{ 8 }$$

The shielding factor ${f}_{{\rm{s}},\mathrm{CO}}({N}_{\mathrm{CO}},{N}_{{{\rm{H}}}_{2}})$ is interpolated from Table 5 in Visser et al. (2009) for a given column density of CO and H₂. This accounts for both the CO self-shielding and the shielding of CO by H₂. In Visser et al. (2009), Table 5, the velocity dispersion of CO is assumed to be $b(\mathrm{CO})=0.3\,\mathrm{km}\,{{\rm{s}}}^{-1}$,⁴ the excitation temperature ${T}_{\mathrm{ex}}(\mathrm{CO})=5{\rm{K}}$, and $N{(}^{12}\mathrm{CO})/N{(}^{13}\mathrm{CO})=69.$

Lastly, the photodissociation rate of C is

$$\begin{eqnarray}&&{R}_{\mathrm{thick},{\rm{C}}}=\chi {R}_{\mathrm{thin},{\rm{C}}}{f}_{\mathrm{dust}}{f}_{{\rm{s}},{\rm{C}}}({N}_{{\rm{C}}},{N}_{{{\rm{H}}}_{2}}).\end{eqnarray} \tag{ 9 }$$

We adopt the treatment in Tielens & Hollenbach (1985). The shielding factor ${f}_{{\rm{s}},{\rm{C}}}({N}_{{\rm{C}}},{N}_{{{\rm{H}}}_{2}})=\exp (-{\tau }_{{\rm{C}}}){f}_{{\rm{s}},{\rm{C}}}({{\rm{H}}}_{2})$, where ${\tau }_{{\rm{C}}}=1.6\times {10}^{-17}{N}_{{\rm{C}}}/{\mathrm{cm}}^{-2}$, ${f}_{{\rm{s}},{\rm{C}}}({{\rm{H}}}_{2})=\exp (-{r}_{{{\rm{H}}}_{2}})/(1+{r}_{{{\rm{H}}}_{2}})$, and ${r}_{{{\rm{H}}}_{2}}=2.8\times {10}^{-22}{N}_{{{\rm{H}}}_{2}}/{\mathrm{cm}}^{-2}$.

The heating and cooling processes included in our models are listed in Table 3, and described in detail in Appendix B. Here we describe the general processes we considered. The time rate of change of the gas thermal energy per H nucleus ${e}_{{\rm{g}},\mathrm{sp}}$ is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}_{{\rm{g}},\mathrm{sp}}}{{dt}}={{\rm{\Gamma }}}_{\mathrm{tot}}-{{\rm{\Lambda }}}_{\mathrm{tot}}.\end{eqnarray} \tag{ 10 }$$

The total heating rate per H nucleus is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{tot}}={{\rm{\Gamma }}}_{\mathrm{PE}}+{{\rm{\Gamma }}}_{\mathrm{CR}}+{{\rm{\Gamma }}}_{{\rm{H}}2\mathrm{gr}}+{{\rm{\Gamma }}}_{{\rm{H}}2\mathrm{pump}}+{{\rm{\Gamma }}}_{{\rm{H}}2\mathrm{diss}},\end{eqnarray} \tag{ 11 }$$

where ${{\rm{\Gamma }}}_{\mathrm{PE}}$ is the heating from the photoelectric effect on dust, ${{\rm{\Gamma }}}_{\mathrm{CR}}$ is the cosmic-ray heating, ${{\rm{\Gamma }}}_{{\rm{H}}2\mathrm{gr}}$ is the heating from ${{\rm{H}}}_{2}$ formation on dust grains, ${{\rm{\Gamma }}}_{{\rm{H}}2\mathrm{pump}}$ is the heating from ${{\rm{H}}}_{2}$ UV pumping, and ${{\rm{\Gamma }}}_{{\rm{H}}2\mathrm{diss}}$ is the heating from ${{\rm{H}}}_{2}$ photodissociation. For typical environments in the cold atomic and molecular ISM, the photoelectric heating dominates at ${A}_{V}\lesssim 1$ and the cosmic-ray heating dominates at ${A}_{V}\gtrsim 1$.

The total cooling rate per H nucleus is

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{tot}}={{\rm{\Lambda }}}_{\mathrm{line}}+{{\rm{\Lambda }}}_{\mathrm{dust}}+{{\rm{\Lambda }}}_{\mathrm{rec}}+{{\rm{\Lambda }}}_{{\rm{H}}2\mathrm{coll}}+{{\rm{\Lambda }}}_{\mathrm{Hion}},\end{eqnarray} \tag{ 12 }$$

where ${{\rm{\Lambda }}}_{\mathrm{line}}$ is the line cooling by atomic and molecular species, ${{\rm{\Lambda }}}_{\mathrm{dust}}$ the cooling from gas−grain collisions, ${{\rm{\Lambda }}}_{\mathrm{rec}}$ the cooling by electron recombination on dust grains, ${{\rm{\Lambda }}}_{{\rm{H}}2\mathrm{coll}}$ the cooling by collisional dissociation of ${{\rm{H}}}_{2}$, and ${{\rm{\Lambda }}}_{\mathrm{Hion}}$ the cooling by collisional ionization of H. The line cooling dominates in typical atomic and diffuse molecular ISM where the gas densities are not high enough for dust cooling to be important ($n\lesssim {10}^{6}\,{\mathrm{cm}}^{-3}$). We include line cooling of O, C, ${{\rm{C}}}^{+}$ fine structure lines, the Lyα line of H, CO rotational lines, and the ${{\rm{H}}}_{2}$ vibration and rotational lines.

The gas temperature is related to the gas energy by (Krumholz 2014)

$$\begin{eqnarray}&&T=\displaystyle \frac{{e}_{{\rm{g}},\mathrm{sp}}}{{c}_{v,{\rm{H}}}}.\end{eqnarray} \tag{ 13 }$$

The specific heat at constant volume ${c}_{v,{\rm{H}}}=\tfrac{1}{2}{k}_{{\rm{B}}}f$, where ${k}_{{\rm{B}}}=1.381\times {10}^{-16}\,\mathrm{erg}\,{{\rm{K}}}^{-1}$ is the Boltzmann constant and $f=\sum {f}_{s}{x}_{s}$ is the degree of freedom in all species per H nucleus. For H, He, He⁺, ${{\rm{H}}}^{+}$, and e, f_s = 3 from translational degrees of freedom. For ${{\rm{H}}}_{2}$, the excitation temperature for rotational and vibrational levels are ${\theta }_{\mathrm{rot}}=170.6\,{\rm{K}}$ and ${\theta }_{\mathrm{vib}}=5984\,{\rm{K}}$ (Tomida et al. 2013). Assuming a fixed ortho-to-para ratio of ${{\rm{H}}}_{2}$, the rotational and vibrational heat capacities of ${{\rm{H}}}_{2}$ are very small at typical molecular cloud temperatures $T\lesssim 50\,{\rm{K}}$. Therefore, we ignore the contribution of rotational and vibrational levels of ${{\rm{H}}}_{2}$ to the specific heat, and simply use

$$\begin{eqnarray}&&{c}_{v,{\rm{H}}}=\displaystyle \frac{3}{2}{k}_{{\rm{B}}}({x}_{{\rm{H}}}+{x}_{{{\rm{H}}}_{2}}+{x}_{{{\rm{H}}}^{+}}+{x}_{\mathrm{He}}+{x}_{{\mathrm{He}}^{+}}+{x}_{{\rm{e}}})\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{3}{2}{k}_{{\rm{B}}}[(1-{x}_{{{\rm{H}}}_{2}})+{x}_{\mathrm{He},\mathrm{tot}}+{x}_{{\rm{e}}}].\end{eqnarray} \tag{ 15 }$$

We consider a simple one-dimensional slab model with uniform density n and FUV radiation incident from one side of the slab, similar to typical 1D PDR models. The incident radiation field is expected to be relatively isotropic over 2π steradians. In Section 3, we use the approach of Wolfire et al. (2010) in order to compare with their PDR code: the isotropic radiation field is approximated by assuming a unidirectional flux incident at an angle of 60° to the normal of the slab surface. Thus, with the incident radiation field strength χ, the field strength at the perpendicular column density N from the slab surface will be ${\chi }_{\mathrm{eff}}=(\chi /2){f}_{\mathrm{shield}}(2N)$, where the factor 1/2 in χ comes from the one-sided slab and the factor 2 in N comes from the angle of 60°. In Section 4, we directly calculate the isotropic radiation field strength by averaging over the incident angels.

In each model, we divide a slab into 10³ logarithmically spaced grid zones in the range of $N={10}^{17}/Z-{10}^{22}/Z\,{\mathrm{cm}}^{-2}$. In each zone, walking inward from the lowest N to the highest N in order, we calculate the equilibrium chemistry and temperature. This is because the chemical states of all exterior zones are needed to calculate the radiation field shielding factor for the next zone.

The evolution of abundances of chemical species ${x}_{s}={n}_{s}/n$ is determined by a set of chemical reactions:

$$\begin{eqnarray}\displaystyle \frac{{{dx}}_{s}}{{dt}} & = & \displaystyle \sum _{i,j}\,(\pm {k}_{2\mathrm{body},{\rm{s}}}^{{ij}}{{nx}}_{i}{x}_{j})+\displaystyle \sum _{i}\,(\pm {k}_{\mathrm{gr},{\rm{s}}}^{i}{{nx}}_{i})\\ & & +\displaystyle \sum _{i}\,(\pm {k}_{\mathrm{cr},{\rm{s}}}^{i}{x}_{i})+\displaystyle \sum _{i}\,(\pm {k}_{\gamma ,{\rm{s}}}^{i}{x}_{i}).\end{eqnarray} \tag{ 16 }$$

Here ${k}_{2\mathrm{body},{\rm{s}}}^{{ij}}$, ${k}_{\mathrm{gr},{\rm{s}}}^{i}$, ${k}_{\mathrm{cr},{\rm{s}}}^{i}$, and ${k}_{\gamma ,{\rm{s}}}^{i}$ are the rates for the two-body, grain surface, cosmic-ray, and photodissociation reactions listed in Tables 1 and 2 that have species s either as a reactant or product. x_i and x_j are the abundance of the corresponding reactant species, and the sign is negative if x_s is a reactant of the chemical reaction and positive if x_s is a product. The chemical reactions in Equation (16) and the energy evolution in Equation (10) give a coupled ODE (ordinary differential equation) system of 13 variables (12 independently calculated species and 1 energy equation). The derived species (H, He, C, O, Si, and e) are not directly calculated from Equation (16) but from the conservation of atomic nuclei and charge. To efficiently solve this set of stiff ODEs, we use the open-source CVODE package (Cohen et al. 1996), which adopts a method of implicit backward differentiation formulas.

Following Krumholz et al. (2011), we investigate how well CO and ${{\rm{H}}}_{2}$ trace cold gas in the ISM with different metallicities. Krumholz et al. (2011) used semi-analytic models to estimate the ${{\rm{H}}}_{2}$ and CO abundances and separately computed the equilibrium temperature for gas cooled either by ${{\rm{C}}}^{+}$ or by CO. Here, we instead use our chemical network and self-consistent cooling to solve for the equilibrium chemical composition and temperature of the gas. We run one-dimensional slab models described in Section 2.3 at a range of densities $n={10}^{1}\mbox{--}{10}^{4}\,{\mathrm{cm}}^{-3}$, and compute the chemical and temperature state of the gas in the A_V–n plane. We also run models with different metallicities, incident radiation field, and cosmic-ray ionization rate to study the dependence of gas properties on these parameters.

It should be noted that in the realistic ISM, the timescales of the cooling and chemical reactions are different, and the dynamical timescales set by turbulence may be shorter than these. Thus, the equilibrium chemical and thermal state does not necessarily apply in the real ISM. This is particularly an issue for ${{\rm{H}}}_{2}$, because the formation time is ${t}_{{{\rm{H}}}_{2}}\approx {10}^{7}\,\mathrm{yr}\,{{\rm{Z}}}_{d}^{-1}{\left(\tfrac{{\rm{n}}}{100{\mathrm{cm}}^{-3}}\right)}^{-1}$ at $T=100\,{\rm{K}}$ (see line 24 in Table 1), whereas typical dynamical timescales in dense gas are (see Equation (36))

$$\begin{eqnarray}&&{t}_{\mathrm{dyn}}=L/{v}_{\mathrm{turb}}{(L)\sim 1\mathrm{Myr}(L/\mathrm{pc})}^{1/2}.\end{eqnarray} \tag{ 17 }$$

For transient clouds produced by turbulent compression, even if the mean shielding column is sufficient for the equilibrium ${{\rm{H}}}_{2}$ abundance to be high, the ${{\rm{H}}}_{2}$ formation timescale may be longer than the cloud lifetime. For example, the criterion ${t}_{{{\rm{H}}}_{2}}\lt {t}_{\mathrm{dyn}}$ is equivalent to

$$\begin{eqnarray}{Z}_{d}^{2}\left(\displaystyle \frac{n}{100\,{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{N}{{10}^{21}\,{\mathrm{cm}}^{-2}}\right) & \gt & 31,\,\mathrm{or}\\ {Z}_{d}{A}_{V}\left(\displaystyle \frac{n}{100\,{\mathrm{cm}}^{-3}}\right) & \gt & 17,\end{eqnarray} \tag{ 18 }$$

which is well above the equilibrium H/${{\rm{H}}}_{2}$ transition in Figure 1.

In high column density self-gravitating clouds with supersonic turbulence, the mass-weighed density $\langle n{\rangle }_{M}$ is proportional to N², $\langle n{\rangle }_{M}\propto {N}^{2}$ (see Equation (50), where the second term dominates over the first term). Combining this with Equation (18) gives

$$\begin{eqnarray}N & \gtrsim & 4.3\times {10}^{21}{Z}_{d}^{-2/3}\,{\mathrm{cm}}^{-2},\,\mathrm{or}\\ {A}_{V} & \gtrsim & 2.3\,{Z}_{d}^{-2/3}\end{eqnarray} \tag{ 19 }$$

for ${t}_{{{\rm{H}}}_{2}}\lt {t}_{\mathrm{dyn}}$. Equation (19) suggests that for a molecular cloud with ${A}_{V}\gtrsim 1$, although molecular hydrogen is difficult to form at an average density $n\sim 100\,{\mathrm{cm}}^{-3}$ (Equation (18)), turbulence can compress most of the gas to reach to a higher density regime and enable the cloud to become molecular in a shorter timescale.

For gravitationally bound clouds that live longer than several flow-crossing times, the dynamical cycling of gas to the unshielded surface may limit the molecular abundance. The timescale for a fluid element in the cloud to travel through the unshielded region is ${t}_{\mathrm{sh}}={L}_{\mathrm{shield}}/{v}_{\mathrm{turb}}({L}_{\mathrm{cloud}})$, where ${L}_{\mathrm{shield}}$ is the length scale for ${{\rm{H}}}_{2}$ to be self-shielded against the FUV radiation. Using Equation (52) in Sternberg et al. (2014) and the ${{\rm{H}}}_{2}$ formation rate on dust (see line 1 in Table 2), the shielding length⁷

$$\begin{eqnarray}&&{L}_{\mathrm{shield}}=\displaystyle \frac{N}{n}\approx 1\,\mathrm{pc}\,\chi {\left(\displaystyle \frac{{\rm{n}}}{100{\mathrm{cm}}^{-3}}\right)}^{-2},\end{eqnarray} \tag{ 20 }$$

where n is the average (volume-weighted) density of the cloud. Adopting ${v}_{\mathrm{turb}}{(L)\sim 1\mathrm{km}{{\rm{s}}}^{-1}(L/\mathrm{pc})}^{1/2}$, this yields:

$$\begin{eqnarray}&&{t}_{\mathrm{sh}}=0.2\,\mathrm{Myr}\,\chi {\left(\displaystyle \frac{n}{100{\mathrm{cm}}^{-3}}\right)}^{-2}{\left(\displaystyle \frac{{L}_{\mathrm{cloud}}}{10\mathrm{pc}}\right)}^{-1/2}.\end{eqnarray} \tag{ 21 }$$

The timescale for ${{\rm{H}}}_{2}$ photodissociation is ${t}_{\mathrm{diss},{{\rm{H}}}_{2}}=1/{k}_{\mathrm{diss},{{\rm{H}}}_{2}}\,=6\times {10}^{2}\,\mathrm{yr}$ (the photodissociation rate ${k}_{\mathrm{diss},{{\rm{H}}}_{2}}=5.6\,\times {10}^{-11}\,{{\rm{s}}}^{-1}$; see line 19 in Table 2). Setting ${t}_{\mathrm{sh}}\lt {t}_{\mathrm{diss},{{\rm{H}}}_{2}}$ gives

$$\begin{eqnarray}&&\displaystyle \frac{1}{\chi }{\left(\displaystyle \frac{n}{100{\mathrm{cm}}^{-3}}\right)}^{2}{\left(\displaystyle \frac{{L}_{\mathrm{cloud}}}{10\mathrm{pc}}\right)}^{1/2}\gt 3\times {10}^{2};\end{eqnarray} \tag{ 22 }$$

when this condition is satisfied, molecules can avoid being destroyed by photodissociation over the timescale to cross the unshielded surface of the cloud. For example, a cloud of size $10\,\mathrm{pc}$ will require its density $n\gt 1.7\times {10}^{3}\,{\mathrm{cm}}^{-3}$ to be fully molecular. Equation (22) shows that even in gravitationally bound clouds, the cycling of gas from the cloud interior to the surface can also limit the ${{\rm{H}}}_{2}$ abundance in low density regimes. If the ${{\rm{H}}}_{2}$ abundance is reduced by exposure to UV, other chemical pathways will also be strongly affected. Therefore, to fully investigate how well CO traces molecular gas, numerical simulations with time-dependent chemistry are needed. Still, we show that our simple models can provide insight into this important but complicated problem.

Figure 5 shows the equilibrium gas temperature in color scale, in comparison to the contours of equilibrium ${{\rm{H}}}_{2}$, CO, and ${{\rm{C}}}^{+}$ abundances. Across a wide range of metallicities $Z=0.2\mbox{--}1$, Figure 5 shows that gas that is mostly molecular is also at $T\lesssim 40\,{\rm{K}}$. However, it is important to note that the coincidence of high ${{\rm{H}}}_{2}$ abundance and low temperature is not causal. For gas in the range of Figure 5, most of the cooling is provided by ${{\rm{C}}}^{+}$, C, and CO. Even if we artificially turn off the formation of ${{\rm{H}}}_{2}$ and CO, the gas can still be cooled to $T\sim 20\,{\rm{K}}$ with C and ${{\rm{C}}}^{+}$ cooling. We discuss this further in Section 4.3. At low metallicity, the minimum density required for gas to become molecular increases. This is because ${{\rm{H}}}_{2}$ forms less efficiently when the dust abundance drops. At low enough n, the equilibrium $x({{\rm{H}}}_{2})$ is low even in very shielded regions. This is due to the destruction of ${{\rm{H}}}_{2}$ by cosmic rays. The gas temperature is similar at ${A}_{V}\lesssim 1$ across $Z=0.2\mbox{--}1$, because both the photoelectric heating and gas cooling are proportional to Z. At ${A}_{V}\gtrsim 1$, and also in low density and low metallicity regions, the cosmic-ray heating dominates, which is independent of metallicity. This leads to an increase of the gas temperature in these regions at low metallicity due to decreased gas cooling.

Zoom In Zoom Out Reset image size

At solar metallicity (Z = 1), Figure 5 shows that the region of high CO abundance (where CO line cooling dominates) traces the lowest temperature $T\lesssim 20\,{\rm{K}}$ gas very well.⁸ At very low metallicities, there is less dust shielding and collisional reactions between metal species occur at a reduced rate, leading to less efficient CO formation, and CO traces the temperature less well.

Figure 6(a) shows that an increase of χ pushes the ${{\rm{H}}}_{2}$-dominated region to higher A_V, and that even ${{\rm{H}}}_{2}$-dominated gas can be at $T\gt 100\,{\rm{K}}$ if χ is large enough. Increasing the cosmic-ray rate requires higher density for the gas to become ${{\rm{H}}}_{2}$ dominated, because cosmic rays can dissociate ${{\rm{H}}}_{2}$ (Figure 6(b)).

Zoom In Zoom Out Reset image size

An increase (decrease) of χ moves the boundary of CO-dominated gas to higher (lower) A_V, because CO only forms in dust-shielded regions (Figure 6(a)). Similarly, an increase (decrease) of ${\xi }_{{\rm{H}}}$ pushes the minimum density for CO up (down), because ${\mathrm{He}}^{+}$ and ${{\rm{C}}}^{+}$ formed in cosmic-ray reactions are harmful for CO formation (Figure 6(b)). However, in both Figures 6(a) and (b), the contours defining the CO-dominated region also clearly delimit the low temperature ($T\lesssim 20\,{\rm{K}}$) gas, except at the highest levels of ${\xi }_{{\rm{H}}}$. The temperature of the bulk of the gas in the CO-dominated regime depends on ${\xi }_{{\rm{H}}}$ but not on χ. This is because when $x(\mathrm{CO})/{x}_{{\rm{C}},\mathrm{tot}}\approx 1$, the heating and cooling are respectively dominated by cosmic-ray ionization and CO rotational transitions. The temperature is obtained by balancing Equation (29) (with Equation (31)) with Equation (33). We suggest that the temperature of CO, especially the optically thin isotopes, can be an indicator to probe the cosmic-ray ionization rate in dense molecular gas. Observations of star-forming regions near the galactic center indicate that the CO gas there is indeed warmer than that in the galactic disk (Güsten et al. 2004; Ao et al. 2013; Bally & Hi-GAL Team 2014). The higher cosmic-ray production rate due to higher SFR in the galactic center may explain the elevated temperature of the CO gas there. However, if only the ${}^{12}\mathrm{CO}\,(J=1-0)$ line is observed, because it is usually optically thick, the luminosity of the line is determined mostly by the temperature at the cloud surface, where photoelectric heating can dominate over cosmic-ray heating. Therefore, as pointed out by Wolfire et al. (1993), the luminosity of the $\mathrm{CO}\,(J=1-0)$ line is mostly determined by the strength of UV radiation, instead of the cosmic-ray ionization rate.

In galactic disks where the warm and cold atomic gas are in pressure equilibrium, the typical density of the cold neutral medium (CNM) is related to the strength of the ambient radiation field approximately by ${n}_{\mathrm{CNM}}\propto \chi$ (Wolfire et al. 2003; Ostriker et al. 2010):

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{CNM}}}{\chi }\approx 23.\end{eqnarray} \tag{ 23 }$$

Krumholz et al. (2009) adopted a similar relation for the density of cold atomic/molecular complexes. Assuming the ambient radiation field strength is proportional to the SFR, then $\chi \propto {{\rm{\Sigma }}}_{\mathrm{SFR}}$, where ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ is the average rate of star formation per unit area in the galactic disk. In equilibrium, ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ is expected to be proportional to the weight of the ISM, which in general depends on both the gas surface density ${{\rm{\Sigma }}}_{\mathrm{gas}}$ and the stellar density (Ostriker et al. 2010; Kim et al. 2011, 2013). In starburst regions where the gas dominates gravity, ${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {{\rm{\Sigma }}}_{\mathrm{gas}}^{2}$ (Ostriker & Shetty 2011). Assuming that ${\xi }_{{\rm{H}}}\propto {{\rm{\Sigma }}}_{\mathrm{SFR}}/{{\rm{\Sigma }}}_{\mathrm{gas}}$ (Ostriker et al. 2010) and ${{\rm{\Sigma }}}_{\mathrm{SFR}}/{{\rm{\Sigma }}}_{\mathrm{gas}}\propto {{\rm{\Sigma }}}_{\mathrm{gas}}\propto \sqrt{\chi }$, we consider the case in which

$$\begin{eqnarray}&&{\xi }_{{\rm{H}}}={10}^{-16}\sqrt{\chi }\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}.\end{eqnarray} \tag{ 24 }$$

Using $\chi \propto n$ and ${\xi }_{{\rm{H}}}\propto \sqrt{\chi }$ based on Equations (23) and (24), Figure 7 plots the temperature of the gas in the n–N plane. Comparing Figures 5 and 7 shows clearly that it is very difficult to form CO in metal-poor gas under typical diffuse-ISM conditions, implying gas densities must be enhanced by turbulence or gravity to be well above typical values in the CNM for CO to be present. However, at Z = 1, most of the carbon is in CO at ${A}_{V}\gtrsim 2$ and $n\gtrsim 100\,{\mathrm{cm}}^{-3}$. Thus, in regions of galaxies where A_V and n are high enough, gas need not be in gravitationally bound clouds to be molecular or to emit strongly in CO (see also Elmegreen 1993). This implies that CO emission in the galactic center and starburst galaxies may arise largely from diffuse gas. Indeed, observations of CO lines in luminous infrared galaxies by Papadopoulos et al. (2012) found that most of the CO emission can come from warm and diffuse gas in these starburst environments. Figure 7 also shows that in metal-poor galaxies, C could be the most important tracer for shielded gas, as also noted by Glover & Clark (2016). However, the gas that is traced by C is not necessarily primarily molecular ${{\rm{H}}}_{2}$.

Zoom In Zoom Out Reset image size

We have also found a simple fit for the contour of $x(\mathrm{CO})/{x}_{{\rm{C}},\mathrm{tot}}=0.5$ in the plane of n – A_V, given values of the incident FUV radiation field, cosmic-ray ionization rate, and metallicity:

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{crit},\mathrm{CO}}}{{\mathrm{cm}}^{-3}}={(4\times {10}^{3}Z{\xi }_{{\rm{H}},16}^{-2})}^{{\chi }_{\mathrm{CO}}^{1/3}}\left(\displaystyle \frac{50{\xi }_{{\rm{H}},16}}{{Z}^{1.4}}\right),\end{eqnarray} \tag{ 25 }$$

where ${n}_{\mathrm{crit},\mathrm{CO}}$ is the critical value above which $x(\mathrm{CO})/{x}_{{\rm{C}},\mathrm{tot}}\,\gt 0.5;$ ${\xi }_{{\rm{H}},16}={\xi }_{{\rm{H}}}/({10}^{-16}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1});$ ${\chi }_{\mathrm{CO}}=\chi \exp (-{\gamma }_{\mathrm{CO}}{A}_{V})$ is the effective radiation field for CO photodissociation accounting for dust attenuation and ${\gamma }_{\mathrm{CO}}=3.53$ is the dust-shielding factor of CO given in Table 2. Figures 5–7 and 9 shows the fit in Equation (25) (white dashed line) against the true contour of $x(\mathrm{CO})/{x}_{{\rm{C}},\mathrm{tot}}\,=0.5$ (magenta dashed line), demonstrating the good agreement between the two. We note that this fit is only tested here to be applicable in the range of $n\approx 10\mbox{--}{10}^{4}\,{\mathrm{cm}}^{-3}$, ${\xi }_{{\rm{H}}}\approx {10}^{-17}\mbox{--}{10}^{-15}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$, and $Z=0.2\mbox{--}1$.

How sensitive is the gas temperature to the detailed chemistry? Glover & Clark (2012b) showed that in simulations of turbulent clouds, as long as the gas is shielded, atomic gas can reach similarly low temperatures as in molecular gas. Glover & Clark (2012a) also found that the gas temperature in dense clouds can be very similar even if the CO and C abundances produced by different chemical networks vary by more than an order of magnitude.

To directly look into this question, we artificially force the gas into different chemical states by setting the formation rates of CO, C, and ${{\rm{H}}}_{2}$ to zero (successively). Figure 8 shows that the equilibrium gas temperatures are very similar with completely different chemical states. C and ${{\rm{C}}}^{+}$ can cool the gas just as well as CO.⁹ ${{\rm{H}}}_{2}$ only affects the temperature slightly by changing the specific heat of the gas. This confirms the conclusions in Glover & Clark (2012b): any correlation between molecular gas and low temperatures is likely to be coincidental instead of causal. Molecular gas traces low temperatures simply because molecules form in dense and shielded regions that are more efficient at cooling and less exposed to heating.

Zoom In Zoom Out Reset image size

Krumholz et al. (2011) argued that ${{\rm{H}}}_{2}$ traces star formation because ${{\rm{H}}}_{2}$ forms in the same regions where the gas is more gravitationally unstable, or, quantitatively speaking, where the critical Bonnor–Ebert mass ${M}_{\mathrm{BE}}=1.2{c}_{s}^{3}{\left(\tfrac{{\pi }^{3}}{{G}^{3}\rho }\right)}^{1/2}=1.2{\left(\tfrac{{\pi }^{3}{k}^{3}{T}^{3}}{{G}^{3}n{\mu }^{4}}\right)}^{1/2}$ (Bonnor 1956; Ebert 1955) becomes low.

In Figure 9, we show the critical Bonnor–Ebert mass as a function of gas density and column similar to Krumholz et al. (2011), with the improvement that chemistry and temperature are calculated self-consistently. As noted by Krumholz et al. (2011), the contours of the 50% and 90% ${{\rm{H}}}_{2}$ abundances are somewhat similar to those of the Bonnor–Ebert mass of the gas. However, we note that the Bonnor–Ebert mass near these contours is ${M}_{\mathrm{BE}}={10}^{2}\mbox{--}{10}^{3}\,{M}_{\odot }$, much larger than the mass of individual stars. In some regions, the Bonnor–Ebert mass is even larger than the available mass in the clouds: for example, at $n=1000\,{\mathrm{cm}}^{-3}$ and A_V = 1, the total mass of the cloud $M\sim {R}^{3}{{nm}}_{{\rm{H}}}\sim {(N/n)}^{3}{{nm}}_{{\rm{H}}}\sim 5\,{M}_{\odot }$, and is much smaller than the local Bonnor–Ebert mass ${M}_{\mathrm{BE}}\sim 100\,{M}_{\odot }$. This suggests that the correlation of ${{\rm{H}}}_{2}$ with star formation is not a sufficient pre-condition, but a fairly non-specific coincidence. To reach Bonnor–Ebert masses comparable to those of individual stars, much higher densities are needed than the densities required for high equilibrium ${{\rm{H}}}_{2}$ abundance.¹⁰

Zoom In Zoom Out Reset image size

The ultimate test for any chemistry model is the comparison with observations. The one-sided equilibrium slab models that we have used to test our network are extremely simple and cannot be expected to agree in detail with the chemical state in the real ISM, which is characterized by complex time-dependent dynamics and morphological structure, both the result of turbulence. Nevertheless, for illustrative purposes it is useful to compare the abundances obtained with our network (for idealized slab geometry) with the observed ISM abundances. To do so, we compiled observations of ${{\rm{H}}}_{2}$, CO, C, ${{\rm{C}}}^{+}$, ${\mathrm{CH}}_{x}$, and ${\mathrm{OH}}_{x}$ abundances in the literature derived from UV and optical absorption spectra along different sight lines: ${{\rm{H}}}_{2}$, C, and ${{\rm{C}}}^{+}$ observations compiled by Wolfire et al. (2008)¹¹ ; ${{\rm{H}}}_{2}$ and CO observations in Rachford et al. (2002) and Sheffer et al. (2008); ${{\rm{H}}}_{2}$ observations in Rachford et al. (2009); ${{\rm{H}}}_{2}$, C, and CO observations in Burgh et al. (2010 hereafter BFJ2010); ${\mathrm{CH}}_{x}$ ¹² abundance in Sheffer et al. (2008) and Crenny & Federman (2004) (hereafter CF2004); and ${\mathrm{OH}}_{x}$ ¹³ abundance in Weselak et al. (2009, 2010). The UV spectroscopy data in these observations mainly come from the Far Ultraviolet Spectroscopic Explorer (FUSE), the Space Telescope Imaging Spectrograph (STIS) on board the Hubble Space Telescope, the UVES spectrograph at the European Southern Observatory (ESO), and the Copernicus survey. The optical spectra data used to derive $\mathrm{CH}$ and CH⁺ abundances in Sheffer et al. (2008) are obtained at the McDonald and European Southern observatories. For most sight lines, the total column N is derived from reddening data using N/E(B − V) = 5.8 × 10²¹ H cm⁻² mag⁻¹ (Bohlin et al. 1978; Rachford et al. 2009), except for sight lines that have H abundances directly observed by Lyα absorption in Bohlin et al. (1978; included in the compilation of observations by Wolfire et al. 2008).

We construct a simple one-sided slab model with thermal and chemical equilibrium, as described in Section 4.1, to compare with the observational data. The fiducial model has $\chi =1$, ${\xi }_{{\rm{H}}}=2\times {10}^{-16}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$, and constant densities of n = 10, 100, and $1000\,{\mathrm{cm}}^{-3}$. We also explore the effects of varying different parameters such as χ, ${\xi }_{{\rm{H}}}$, the efficiency of grain-assisted recombination, and gas-phase carbon abundance.

Figure 10 shows the comparison between the observations and our slab cloud model in column densities of ${{\rm{H}}}_{2}$, CO, C, ${{\rm{C}}}^{+}$ and ${\mathrm{CH}}_{x}$, and ${\mathrm{OH}}_{x}$. For the H to ${{\rm{H}}}_{2}$ transition in the left panel, low density $n\sim 10\,{\mathrm{cm}}^{-2}$ is required to match equilibrium ${{\rm{H}}}_{2}$ abundances with observations. This is consistent with the analysis by Bialy et al. (2015, 2017) on the H to ${{\rm{H}}}_{2}$ transition layers in the Perseus molecular cloud and star-forming region W43.¹⁴

Zoom In Zoom Out Reset image size

Comparing the right panel with the left, there is much less dispersion in the observed C, CO, and ${\mathrm{CH}}_{x}$ abundances when plotted against the ${{\rm{H}}}_{2}$ column density ${N}_{{{\rm{H}}}_{2}}$ instead of the total column N or A_V. This is likely due to foreground or background contaminations: low density atomic ISM is much more "diffuse" in spatial distribution than the relatively "clumpy" molecular ISM. Therefore, the low density atomic clouds may contribute a significant fraction of the total column/visual extinction in a given sight line, without contributing much to the total column of chemical species such as ${{\rm{H}}}_{2}$, C, CO, ${\mathrm{CH}}_{x}$, and ${\mathrm{OH}}_{x}$ which mainly form in higher density regions.

In spite of the extreme idealizations adopted, there is generally a good agreement between our model and the observations in Figure 10, especially in the right panel (when plotted against ${N}_{{{\rm{H}}}_{2}}$): we can successfully reproduce the range of observed abundances of CO, ${\mathrm{CH}}_{x}$, and ${\mathrm{OH}}_{x}$ with densities between $100\,{\mathrm{cm}}^{-3}$ and $1000\,{\mathrm{cm}}^{-3}$. It is especially remarkable that CO and ${\mathrm{CH}}_{x}$ abundances agree over one to two orders of magnitude of ${N}_{{{\rm{H}}}_{2}}$.

However, there is one significant discrepancy between the simple model and observations: the predicted C abundance is too high at a given A_V or ${N}_{{{\rm{H}}}_{2}}$ by about an order of magnitude (see Figure 10). This is more clearly shown in the left panel of Figure 11 when plotted with the ratio of CO/C (see also Figure 4 in BFJ2010). To investigate the dependence of CO/C on model parameters, we vary the incident radiation field strength χ, cosmic-ray ionization rate ${\xi }_{{\rm{H}}}$, the efficiency of grain-assisted recombination of ions, and the gas-phase carbon abundance. The summary of our results is shown in the left panel of Figure 11. Variations in χ and ${\xi }_{{\rm{H}}}$ tend to change the C and CO abundances in the same direction, giving a very similar CO/C ratio to that of the fiducial model. Excluding the grain-assisted recombination of ions does not solve the problem either: it does lower the C abundance at a given A_V, but at the same time, the CO abundance is much lower as well, and the CO/C ratio is still too low (see also the left panel of Figure 26).

Zoom In Zoom Out Reset image size

Note that our conclusion is different from the conclusion of in Liszt (2011), which stated that the absence of grain-assisted recombination can reproduce the high CO/C observed in BFJ2010. We believe this is because Liszt (2011) assumed a constant ${x}_{{\mathrm{HCO}}^{+}}\sim 3\times {10}^{-9}$ in their chemistry model. Without grain-assisted recombination, the HCO⁺ abundance is much lower at ${x}_{{\mathrm{HCO}}^{+}}\lesssim {10}^{-10}$, leading to much less efficient CO formation than calculated in Liszt (2011).¹⁵

The only variation of our models that can reproduce the observed CO/C ratio is the one with gas-phase carbon abundance depleted by an order of magnitude relative to the fiducial model, i.e., ${{\rm{C}}}_{\mathrm{tot}}/{\rm{H}}=0.1{({{\rm{C}}}_{\mathrm{tot}}/{\rm{H}})}_{\mathrm{fiducial}}=1.6\times {10}^{-5}$. This is also found in BFJ2010 (their Figure 4). However, although carbon depletion can be caused by the formation of grains such as polycyclic aromatic hydrocarbons (PAHs), a depletion factor as extreme as 0.1 is unlikely. Even the most carbon-depleted sight lines observed in Sofia et al. (2011) and Parvathi et al. (2012) have ${{\rm{C}}}_{\mathrm{tot}}/{\rm{H}}\gt 6\times {10}^{-5}$ (derived from ${{\rm{C}}}^{+}$ absorption measurements), and many of these sight lines are the same as in Burgh et al. (2010). Moreover, as shown in the right panel of Figure 11, the ratio of CO/${\mathrm{CH}}_{x}$ in the model with carbon depletion is an order of magnitude higher than the observed values, whereas other models reproduce CO/${\mathrm{CH}}_{x}$ consistent with observations. Therefore, we conclude that the depletion of gas-phase carbon is unlikely to be the solution for this problem.

Overpredicting C abundance is also an issue among other PDR models. For example, BFJ2010 compared their observations with the chemistry model by van Dishoeck & Black (1988) and found very similar results. Glover et al. (2010) also predicted $\mathrm{CO}/{\rm{C}}\ll 1$ for ${A}_{V}\lt 2$ clouds with their extensive chemical network (see their Figures 1 and 2). Bensch et al. (2003) observed C and CO emission in the translucent cloud MCLD 123.5+24.9, and also found $\mathrm{CO}/{\rm{C}}\gtrsim 1$. They used the PDR code by Stoerzer et al. (1996) to model their observations, and found that their models also tend to produce a CO/C ratio that is too low, unless they assume a cloud structure of high density ($n\sim {10}^{4}\,{\mathrm{cm}}^{-3}$) clumps embedded in a low density medium.

It is worth emphasizing that observations compiled here are all for ${A}_{V}\lesssim 1$ diffuse and translucent clouds, for which UV or optical absorption observations may be used. For these clouds, most carbon (>90%) is still in the form of ${{\rm{C}}}^{+}$. The physical properties of these clouds are very different from GMCs and smaller dark molecular clouds, which typically has ${A}_{V}\gg 1$ and most carbon is in the form of CO.

There are some potential solutions to the mismatch of C abundance between equilibrium PDR models and observations. (1) Non-equilibrium chemistry. As previously noted, the photoionization and photodissociation timescales are generally short compared to collisional reactions or grain-assisted reactions (see Tables 1 and 2). Dynamical effects such as limited cloud lifetime and turbulent cycling of gas from a cloud's interior to its surface can lower the abundances of species formed in low-shielding regions such as C. Furthermore, C tends to form in less dense and shielded regions than CO, and can be more subject to these dynamical effects, which may bring its abundance far from equilibrium. (2) Improvement of chemical networks and reaction rates. It is possible that the chemical pathways to form C and CO in translucent clouds are not well understood, causing discrepancy between chemical models and observations. For example, the composition of dust grains and the rates for grain surface reactions are still uncertain. In addition, turbulent dissipation (Godard et al. 2014) might convert the dominant reservoir of carbon (in ${{\rm{C}}}^{+}$) to CO without much production of C. In the future, 3D realistic ISM simulations with time-dependent chemistry will give us more insight into the role of non-equilibrium chemistry.

We have added or modified the following reactions to the NL99 network:

• 1.  
  ${{\rm{H}}}_{3}^{+}$ destruction channel

  $$\begin{eqnarray*}&&{{\rm{H}}}_{3}^{+}\ +{\rm{e}}\ \to \ 3{\rm{H}}\end{eqnarray*}$$

  in addition to

  $$\begin{eqnarray*}&&{{\rm{H}}}_{3}^{+}\ +{\rm{e}}\ \to {{\rm{H}}}_{2}\ +{\rm{H}}\end{eqnarray*}$$

  in the original NL99 network, to be more consistent with the hydrogen network adopted in Glover et al. (2010) and Indriolo et al. (2007).
• 2.  
  Grain-assisted recombination of ${{\rm{C}}}^{+}$ and He⁺:

  $$\begin{eqnarray*}{{\rm{C}}}^{+}+{\rm{e}}+\mathrm{gr} & \to & {\rm{C}}+\mathrm{gr}\\ {\mathrm{He}}^{+}+{\rm{e}}+\mathrm{gr} & \to & \mathrm{He}+\mathrm{gr}.\end{eqnarray*}$$

  These two grain surface reactions are the main channels for ${{\rm{C}}}^{+}$ and He⁺ recombination at solar metallicity. Glover & Clark (2012a) considered these two reactions, and concluded that they made very little difference in determining where CO would form. However, we found that using the updated cosmic-ray ionization rate ${\xi }_{{\rm{H}}}=2\times {10}^{-16}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$ instead of their rate ${\xi }_{{\rm{H}}}\,={10}^{-17}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$, these reactions, especially the recombination of ${{\rm{C}}}^{+}$ on dust grains, are essential for CO formation at moderate densities $n\lesssim 1000\,{\mathrm{cm}}^{-3}$. This is because electrons from cosmic-ray ionization of H and C inhibit CO formation, while He⁺ formed by cosmic rays destroys CO.
• 3.  
  Cosmic-ray ionization of H and ${{\rm{H}}}_{2}$ by secondary electrons. The measurement of the cosmic-ray ionization rate ${\xi }_{{\rm{H}}}=2.0\times {10}^{-16}\,{{\rm{s}}}^{-1}\,{{\rm{H}}}^{-1}$ in Indriolo et al. (2007, ${\xi }_{p}$ in their paper) only accounts for the initial (primary) cosmic-ray ionization. The secondary electrons created by the first ionization event can again ionize H and ${{\rm{H}}}_{2}$. The secondary ionization rate depends on factors such as electron abundance (Dalgarno et al. 1999). Here we apply a simple approach motivated by Glassgold & Langer (1974), which has the total ionization rate of H 1.5 times the primary rate in atomic regions and 1.15 times in molecular regions. In regions mixed with atomic and molecular gas, we simply scale the total rate with H and ${{\rm{H}}}_{2}$ abundance as shown in reactions 6 and 7 in Table 2. The cosmic-ray ionization rate of ${{\rm{H}}}_{2}$ is simply twice that of H.
• 4.  
  Removed photodissociation of HCO⁺. We find that the destruction of HCO⁺ is dominated by the reaction ${\mathrm{HCO}}^{+}+{\rm{e}}\to \mathrm{CO}+{\rm{H}}$, and the photodissociation of HCO⁺ has very little effect on the whole chemistry network.
• 5.  
  He⁺ destruction in the ${{\rm{H}}}_{2}$ region:

  $$\begin{eqnarray*}{\mathrm{He}}^{+}+{{\rm{H}}}_{2} & \to & {{\rm{H}}}_{2}^{+}+\mathrm{He}\\ {\mathrm{He}}^{+}+{{\rm{H}}}_{2} & \to & {{\rm{H}}}^{+}+\mathrm{He}+{\rm{H}}.\end{eqnarray*}$$

  The rates of these two reactions are similar to the rate of direct recombination of He⁺ with electrons, and can be higher than the grain surface recombination of He⁺ abundance at low metallicities.
• 6.  
  CH destruction:

  $$\begin{eqnarray*}&&\mathrm{CH}+{\rm{H}}\to {{\rm{H}}}_{2}+{\rm{C}}.\end{eqnarray*}$$

  This is one of the main destruction channels of CH in shielded regions.
• 7.  
  Cosmic-ray ionization and cosmic-ray- induced photodissociation of ${\rm{C}}$ and CO:

  $$\begin{eqnarray*}\mathrm{cr}+{\rm{C}} & \to & {{\rm{C}}}^{+}+{\rm{e}}\\ \mathrm{cr}+\mathrm{CO} & \to & {\mathrm{CO}}^{+}+{\rm{e}}\\ {\gamma }_{\mathrm{cr}}+{\rm{C}} & \to & {{\rm{C}}}^{+}+{\rm{e}}\\ {\gamma }_{\mathrm{cr}}+\mathrm{CO} & \to & {\rm{C}}+{\rm{O}}.\end{eqnarray*}$$

  For $\mathrm{cr}+\mathrm{CO}\to {\mathrm{CO}}^{+}+{\rm{e}}$, we assume that ${\mathrm{CO}}^{+}$ reacts quickly with H/${{\rm{H}}}_{2}$ to form HCO⁺ and implement the reaction as $\mathrm{cr}+\mathrm{CO}+{\rm{H}}\to {\mathrm{HCO}}^{+}+{\rm{e}}$ with the same rate. Similar to Clark & Glover (2015), the rates for these reactions are assumed to be proportional to the cosmic-ray ionization rate of atomic hydrogen ${\xi }_{{\rm{H}}}$. The scaling factors for cosmic-ray ionization ${\xi }_{{\rm{C}},\mathrm{CO}}/{\xi }_{{\rm{H}}}$ are calculated with the rates in McElroy et al. (2013), and the scaling factors for cosmic-ray-induced photoionization are adopted from Gredel et al. (1987). We also scale the cosmic-ray-induced photoionization to the ${{\rm{H}}}_{2}$ fraction $2n({{\rm{H}}}_{2})/(2n({{\rm{H}}}_{2})+2n({\rm{H}}))$, because the ionizing photons come from the UV light induced by the primary ionization and excitation of ${{\rm{H}}}_{2}$, and would not be present in regions with only atomic hydrogen. We have also tested the inclusion of cosmic-ray-induced photoreactions of CH and OH using the rates in Gredel et al. (1989), but found they made very little difference for CO formation.
• 8.  
  ${\mathrm{OH}}_{x}+{\rm{O}}\to 2{\rm{O}}+{\rm{H}}$. This reaction can be important for ${\mathrm{OH}}_{x}$ destruction in dense and shielded regions, where photodissociation of ${\mathrm{OH}}_{x}$ is less efficient. See ${\mathrm{OH}}_{x}$ pseudo-reactions in Appendix A.3.
• 9.  
  ${\mathrm{He}}^{+}+{\mathrm{OH}}_{x}\to {{\rm{O}}}^{+}+\mathrm{He}+{\rm{H}}$: This can be the major channel for He⁺ destruction in shielded regions. See Appendix A.3.
• 10.  
  ${{\rm{C}}}^{+}+{{\rm{H}}}_{2}+{\rm{e}}\to {\rm{C}}+2{\rm{H}}$. See the ${{\rm{C}}}^{+}+{{\rm{H}}}_{2}$ reaction in Appendix A.2. Note that this is not a three-body reaction but a pseudo-reaction representing a two-body reaction followed by a recombination.
• 11.  
  ${{\rm{H}}}_{2}{{\rm{O}}}^{+}$ destruction by reacting with electrons:

  $$\begin{eqnarray*}&&{{\rm{H}}}_{2}{{\rm{O}}}^{+}+{\rm{e}}\to 2{\rm{H}}+{\rm{O}}.\end{eqnarray*}$$

  In regions where electron abundance is relatively high, this can be limiting for the abundance of ${\mathrm{OH}}_{x}$ species. This is applied implicitly as the branching ratio for ${\mathrm{OH}}_{x}$ formation in reactions 2–5 in Table 1. See Appendix A.3 for details.
• 12.  
  ${{\rm{O}}}^{+}$ species and reactions:

  $$\begin{eqnarray*}{{\rm{H}}}^{+}+{\rm{O}} & \to & {{\rm{O}}}^{+}+{\rm{H}}\\ {{\rm{O}}}^{+}+{\rm{H}} & \to & {{\rm{H}}}^{+}+{\rm{O}}.\end{eqnarray*}$$

  In high temperature regions, the first reaction has a higher rate than the second one (its reverse reaction). ${{\rm{O}}}^{+}$ formed by change exchange with ${{\rm{H}}}^{+}$ can react with ${{\rm{H}}}_{2}$, and subsequently from OH⁺. See Appendix A.3.
• 13.  
  $\mathrm{Si}$, ${\mathrm{Si}}^{+}$ species and reactions:

  $$\begin{eqnarray*}{\mathrm{Si}}^{+}+{\rm{e}} & \to & \mathrm{Si}\\ {\mathrm{Si}}^{+}+{\rm{e}}+\mathrm{gr} & \to & \mathrm{Si}+\mathrm{gr}\\ {\gamma }_{\mathrm{cr}}+\mathrm{Si} & \to & {\mathrm{Si}}^{+}+{\rm{e}}\\ \gamma +\mathrm{Si} & \to & {\mathrm{Si}}^{+}+{\rm{e}}.\end{eqnarray*}$$

  The NL99 network includes the species ${\rm{M}}$ to represent the metals besides ${\rm{C}}$ and O. However, the metal abundance they are using, ${x}_{{\rm{M}},\mathrm{tot}}=2\times {10}^{-7}$, is much lower than the observed abundances in the solar neighborhood: the most abundant gas-phase metals are ${\rm{S}}$ and $\mathrm{Si}$, with ${x}_{{\rm{S}},\mathrm{tot}}=3.5\times {10}^{-6}$ (Jenkins 2009) and ${x}_{\mathrm{Si},\mathrm{tot}}=1.7\times {10}^{-6}$ (Cardelli et al. 1994). Because the dust-shielding factor ${\gamma }_{\mathrm{Si}}\lt {\gamma }_{{\rm{C}}}$ (see Equation (2) and Table 2), at ${A}_{V}\sim 1$, ${\mathrm{Si}}^{+}$ can be the main agent for providing electrons where C is already shielded from photoionization (see, e.g., Figure 14). Although the total abundance of ${\rm{S}}$ is higher than Si, we find that including S and ${{\rm{S}}}^{+}$ has very little effect on the electron or CO abundances, since ${\gamma }_{{\rm{S}}}\approx {\gamma }_{{\rm{C}}}$. Therefore, we only include $\mathrm{Si}$ and ${\mathrm{Si}}^{+}$ for the metals.

The pseudo-species ${\mathrm{CH}}_{x}$ includes CH, CH₂, CH⁺, ${\mathrm{CH}}_{2}^{+}$, and ${\mathrm{CH}}_{3}^{+}$. The creation of ${\mathrm{CH}}_{x}$ is from two reactions. The first one is

$$\begin{eqnarray*}&&{{\rm{H}}}_{3}^{+}+{\rm{C}}\to {\mathrm{CH}}_{x}+{{\rm{H}}}_{2}.\end{eqnarray*}$$

This represents the sum of two reactions: ${{\rm{H}}}_{3}^{+}+{\rm{C}}\,\to {\mathrm{CH}}^{+}+{{\rm{H}}}_{2}$ and ${{\rm{H}}}_{3}^{+}+{\rm{C}}\to {\mathrm{CH}}_{2}^{+}+{\rm{H}}$. We use the sum of both rates for this pseudo-reaction. CH⁺ reacts with ${{\rm{H}}}_{2}$ to form ${\mathrm{CH}}_{2}^{+}$ and H, which again react with ${{\rm{H}}}_{2}$ to form ${\mathrm{CH}}_{3}^{+}$ and H. The second one,

$$\begin{eqnarray*}&&{{\rm{C}}}^{+}+{{\rm{H}}}_{2}\to {\mathrm{CH}}_{x}+{\rm{H}},\end{eqnarray*}$$

represents the reaction of ${{\rm{C}}}^{+}+{{\rm{H}}}_{2}\to {\mathrm{CH}}_{2}^{+}$. ${\mathrm{CH}}_{2}^{+}$ quickly reacts with ${{\rm{H}}}_{2}$ and turns into ${\mathrm{CH}}_{3}^{+}$ and H, which then reacts with an electron. 70% of the reaction ${\mathrm{CH}}_{3}^{+}+{\rm{e}}$ gives CH or CH₂, and 30% turns back into ${\rm{C}}$. Therefore, we use the rate in Wakelam et al. (2010) multiplied by 0.7 for ${{\rm{C}}}^{+}+{{\rm{H}}}_{2}\to {\mathrm{CH}}_{x}+{\rm{H}}$, and by 0.3 for ${{\rm{C}}}^{+}+{{\rm{H}}}_{2}+{\rm{e}}\to {\rm{C}}+2{\rm{H}}$.

The recombinations of CH⁺, ${\mathrm{CH}}_{2}^{+}$, and ${\mathrm{CH}}_{3}^{+}$ with electrons form CH and CH₂, which can be destroyed in three pathways:

$$\begin{eqnarray*}{\mathrm{CH}}_{x}+{\rm{O}} & \to & \mathrm{CO}+{\rm{H}}\\ {\mathrm{CH}}_{x}+{\rm{H}} & \to & {\rm{C}}+{{\rm{H}}}_{2}\\ \gamma +{\mathrm{CH}}_{x} & \to & {\rm{C}}+{\rm{H}}.\end{eqnarray*}$$

For the reaction ${\mathrm{CH}}_{x}+{\rm{O}}\to \mathrm{CO}+{\rm{H}}$, CH and CH₂ have similar rates. Here we use the rate of CH. For the reactions ${\mathrm{CH}}_{x}+{\rm{H}}\to {\rm{C}}+{{\rm{H}}}_{2}$ and $\gamma +{\mathrm{CH}}_{x}\to {\rm{C}}+{\rm{H}}$, CH₂ will react with H to form CH and with a photon to form CH and ${\mathrm{CH}}_{2}^{+}$. Therefore, we again use the rate assuming ${\mathrm{CH}}_{x}$ are all in CH.

The pseudo-species ${\mathrm{OH}}_{x}$ represents OH, H₂O, OH⁺, H₂O⁺, and H₃O⁺. The formation of ${\mathrm{OH}}_{x}$ is mainly through

$$\begin{eqnarray}&&{{\rm{H}}}_{3}^{+}+{\rm{O}}\to {\mathrm{OH}}_{x}+{{\rm{H}}}_{2}.\end{eqnarray} \tag{ 26 }$$

This represents two reactions, the formation of OH⁺ and H₂O⁺. The sum of both rates gives the rate of this pseudo-reaction. Another channel for ${\mathrm{OH}}_{x}$ formation is

$$\begin{eqnarray}&&{{\rm{O}}}^{+}+{{\rm{H}}}_{2}\to {\mathrm{OH}}_{x}+{\rm{H}},\end{eqnarray} \tag{ 27 }$$

and this represents the reaction ${{\rm{O}}}^{+}+{{\rm{H}}}_{2}\to {\mathrm{OH}}^{+}+{\rm{H}}$. ${{\rm{O}}}^{+}$ mainly comes from the charge exchange ${\rm{O}}+{{\rm{H}}}^{+}\to {{\rm{O}}}^{+}+{\rm{H}}$.

OH⁺ formed by the two channels above reacts with ${{\rm{H}}}_{2}$ to form H₂O⁺. The H₂O⁺ molecule has two fates: it can again react with ${{\rm{H}}}_{2}$ to form H₃O⁺:

$$\begin{eqnarray*}&&{{\rm{H}}}_{2}{{\rm{O}}}^{+}+{\rm{H}}2\to {{\rm{H}}}_{3}{{\rm{O}}}^{+}+{\rm{H}},\,{k}_{{{\rm{H}}}_{2}{{\rm{O}}}^{+},{{\rm{H}}}_{2}}=6.0\times {10}^{-10}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1},\end{eqnarray*}$$

and then H₃O⁺ combines with electrons, eventually forming more stable OH and H₂O. Or, the H₂O⁺ molecule can be destroyed by electrons:

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{{\rm{O}}}^{+}+{\rm{e}}\to 2{\rm{H}}+{\rm{O}},\,{k}_{{{\rm{H}}}_{2}{{\rm{O}}}^{+},{\rm{e}}}=5.3\times {10}^{-6}{T}^{-0.5}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 28 }$$

At low density and low A_V regions where ${{\rm{H}}}_{2}$ abundance is low and electron abundance is high, H₂O⁺ destruction by Equation (28) can limit the formation of ${\mathrm{OH}}_{x}$ by Equations (26) and (27). Therefore, we multiply the reaction rates of Equations (26) and (27) by a branching factor $r={k}_{{{\rm{H}}}_{2}{{\rm{O}}}^{+},{{\rm{H}}}_{2}}{x}_{{{\rm{H}}}_{2}}/({k}_{{{\rm{H}}}_{2}{{\rm{O}}}^{+},{{\rm{H}}}_{2}}{x}_{{{\rm{H}}}_{2}}+{k}_{{{\rm{H}}}_{2}{{\rm{O}}}^{+},{\rm{e}}}{x}_{{\rm{e}}})$, and also add reactions 3 and 5 in Table 1, which is equivalent to the combined reactions of Equations (26) and (27) with Equation (28).

We assume most ${\mathrm{OH}}_{x}$ are in OH, and use the rate of OH reactions for the following pseudo-reactions of ${\mathrm{OH}}_{x}$ destruction:

$$\begin{eqnarray*}{{\rm{C}}}^{+}+{\mathrm{OH}}_{x} & \to & {\mathrm{HCO}}^{+}\\ {\mathrm{OH}}_{x}+{\rm{C}} & \to & \mathrm{CO}+{\rm{H}}\\ \gamma +{\mathrm{OH}}_{x} & \to & {\rm{O}}+{\rm{H}}\\ {\mathrm{OH}}_{x}+{\rm{O}} & \to & 2{\rm{O}}+{\rm{H}}\\ {\mathrm{He}}^{+}+{\mathrm{OH}}_{x} & \to & {{\rm{O}}}^{+}+\mathrm{He}+{\rm{H}}.\end{eqnarray*}$$

For the reaction ${{\rm{C}}}^{+}+{\mathrm{OH}}_{x}\to {\mathrm{HCO}}^{+}$, we use the rate of ${{\rm{C}}}^{+}+\mathrm{OH}\to {\mathrm{CO}}^{+}+{\rm{H}}$, because ${\mathrm{CO}}^{+}$ will quickly react with ${{\rm{H}}}_{2}$ to form HCO⁺. We use the rate of reaction $\mathrm{OH}+{\rm{O}}\to {{\rm{O}}}_{2}+{\rm{H}}$ for ${\mathrm{OH}}_{x}+{\rm{O}}\to 2{\rm{O}}+{\rm{H}}$, assuming most oxygen will be in the form of atomic O.

B.1.1. Cosmic-Ray Ionization

When the gas is ionized by cosmic rays, the primary and secondary electrons thermalize with and heat up the gas. The heating rate per H,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{cr}}=\xi {q}_{\mathrm{cr}},\end{eqnarray} \tag{ 29 }$$

where $\xi ={\xi }_{{\rm{H}}}x({\rm{H}})+{\xi }_{{{\rm{H}}}_{2}}x({{\rm{H}}}_{2})+{\xi }_{\mathrm{He}}x(\mathrm{He})$ and ${q}_{\mathrm{cr}}$ is the energy added to the gas per primary ionization. We use the result fitted by Draine (2011) to the data of Dalgarno & McCray (1972) in atomic regions:

$$\begin{eqnarray}&&{q}_{\mathrm{cr},{\rm{H}}}=6.5+26.4{\left(\displaystyle \frac{x({\rm{e}})}{x(e)+0.07}\right)}^{0.5}\,\mathrm{eV},\end{eqnarray} \tag{ 30 }$$

and the fit by Krumholz (2014) to the data of Glassgold et al. (2012) in molecular regions:

$$\begin{eqnarray}\displaystyle \frac{{q}_{\mathrm{cr},{{\rm{H}}}_{2}}}{\mathrm{eV}}=\left\{\begin{array}{ll}10,\, & \mathrm{log}n\lt 2\\ 10+\displaystyle \frac{3(\mathrm{log}n-2)}{2},\, & 2\leqslant \mathrm{log}n\lt 4\\ 13+\displaystyle \frac{4(\mathrm{log}n-4)}{3},\, & 4\leqslant \mathrm{log}n\lt 7\\ 17+\displaystyle \frac{\mathrm{log}n-7}{3},\, & 7\leqslant \mathrm{log}n\lt 10\\ 18,\, & \mathrm{log}n\geqslant 10\end{array}\right..\end{eqnarray} \tag{ 31 }$$

Similar to Krumholz (2014), we simply assume that the total heating rate can be calculated by summing the heating rates in the atomic and molecular regions weighted by their number fraction:

$$\begin{eqnarray}&&{q}_{\mathrm{cr}}=x({\rm{H}}){q}_{\mathrm{cr},{\rm{H}}}+2x({{\rm{H}}}_{2}){q}_{\mathrm{cr},{{\rm{H}}}_{2}}.\end{eqnarray} \tag{ 32 }$$

B.1.2. Photoelectric Effect on Dust Grains

B.1.3. H₂ Formation on Dust Grains, and UV Pumping of H₂

B.1.4. H₂ Photodissociation

B.2.1. Atomic Lines: ${{\rm{C}}}^{+}$, ${\rm{C}}$, and O Fine Structure Lines, and the Ly$\alpha$ Line

B.2.2. CO Rotational Lines

The CO rotational lines are often optically thick in regions where CO cooling is important. To calculate the cooling rate by CO rotational lines, we use the cooling functions in Omukai et al. (2010, their Appendix C and Table B2). They use the large velocity gradient (LVG) approximation, similar to the approach by Neufeld & Kaufman (1993). The cooling rate per H is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CO}}={L}_{\mathrm{CO}}x(\mathrm{CO}){n}_{{{\rm{H}}}_{2}},\end{eqnarray} \tag{ 33 }$$

where ${L}_{\mathrm{CO}}$ is given in terms of ${n}_{{{\rm{H}}}_{2}}$, $\tilde{N}(\mathrm{CO})$, and T by Equation (B1) in Omukai et al. (2010). In Neufeld & Kaufman (1993) and Omukai et al. (2010), they assume ${{\rm{H}}}_{2}$ is the only species responsible for the collisional excitation of the CO rotational levels, which is appropriate in fully molecular regions. To take into account the collisional excitation by atomic hydrogen and electrons, we replace ${n}_{{{\rm{H}}}_{2}}$ in Equation (33) with an effective number density, following Yan (1997), Meijerink & Spaans (2005), and Glover et al. (2010):

$$\begin{eqnarray}&&{n}_{\mathrm{eff},\mathrm{CO}}={n}_{{{\rm{H}}}_{2}}+\sqrt{2}\left(\displaystyle \frac{{\sigma }_{{\rm{H}}}}{{\sigma }_{{{\rm{H}}}_{2}}}\right)n+\left(\displaystyle \frac{1.3\times {10}^{-8}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}}{{\sigma }_{{{\rm{H}}}_{2}}{v}_{{\rm{e}}}}\right){n}_{{\rm{e}}},\end{eqnarray} \tag{ 34 }$$

where ${\sigma }_{{\rm{H}}}=2.3\times {10}^{-15}\,{\mathrm{cm}}^{2}$, ${\sigma }_{{{\rm{H}}}_{2}}=3.3\times {10}^{-16}{(T/1000{\rm{K}})}^{-1/4}\,{\mathrm{cm}}^{2}$ and ${v}_{{\rm{e}}}=1.03\times {10}^{4}{(T/1{\rm{K}})}^{1/2}\,\mathrm{cm}\,{{\rm{s}}}^{-1}$.

In LVG approximation, the escape probability of a photon emitted by CO is related to the effective column density parameter

$$\begin{eqnarray}&&\tilde{N}(\mathrm{CO})=\displaystyle \frac{n(\mathrm{CO})}{| {\rm{\nabla }}v| },\end{eqnarray} \tag{ 35 }$$

where $n(\mathrm{CO})$ is the local CO number density and ${\rm{\nabla }}v$ is the local gas velocity gradient. In Omukai et al. (2010), they used the approximation $\tilde{N}(\mathrm{CO})=N(\mathrm{CO})/\sqrt{2{kT}/m(\mathrm{CO})}$, where T is the gas temperature, $N(\mathrm{CO})$ the total CO column density, and $m(\mathrm{CO})$ the mass of a single CO molecule. This assumes thermal motions dominate the velocity of CO molecules. However, in realistic molecular clouds, the gas velocity will be largely determined by the supersonic turbulence. For Milky Way molecular clouds, the observed turbulence spectrum obeys the line-width size relation (Larson 1981; Solomon et al. 1987; Heyer & Brunt 2004):

$$\begin{eqnarray}&&v(L)\sim 1\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\left(\displaystyle \frac{L}{\mathrm{pc}}\right)}^{1/2}.\end{eqnarray} \tag{ 36 }$$

In our slab models, we use the velocity determined using Equation (36) to calculate the effective CO column density parameter $\tilde{N}(\mathrm{CO})=N(\mathrm{CO})/v(L)$, where $N(\mathrm{CO})$ is the column density of CO from the local position to the edge of the one-sided slab. This assumes that for thicker clouds, the global turbulent velocity will be larger as implied by Equation (36). We have also tried using a constant velocity $v=1\,\mathrm{km}\,{{\rm{s}}}^{-1}$ to calculate $\tilde{N}$, and the results for the equilibrium temperature in CO-dominated regions are very similar.

Care needs to be taken when the temperature $T\lt 10\,{\rm{K}}$. Omukai et al. (2010) only gives fitting parameters down to $T=10\,{\rm{K}}$, and naive extrapolation may give a finite cooling rate even when the temperature is close to zero. Moreover, CO can freeze-out onto dust grains in regions with high density and $T\lesssim 10\,{\rm{K}}$, leading to depleted CO abundance (Bergin et al. 1995; Acharyya et al. 2007). For simplicity, we artificially turn off CO cooling at temperatures below $10\,{\rm{K}}$. Omukai et al. (2010) also only gives fitting parameters up to $T=2000\,{\rm{K}}$. For $T\gt 2000\,{\rm{K}}$, we simply use the CO cooling rate at $T=2000\,{\rm{K}}$. In reality, CO will likely be destroyed by collisional dissociation and ionization at such high temperatures.

B.2.3. H₂ Rotational and Vibrational Lines

We use the results in Glover & Abel (2008) to calculate cooling by ${{\rm{H}}}_{2}$ rotational and vibrational lines. The cooling rate per ${{\rm{H}}}_{2}$ molecule ${{\rm{\Lambda }}}_{{{\rm{H}}}_{2}}$ is approximated by

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{{\rm{H}}}_{2}}=\displaystyle \frac{{{\rm{\Lambda }}}_{{{\rm{H}}}_{2},\mathrm{LTE}}}{1+{{\rm{\Lambda }}}_{{{\rm{H}}}_{2},\mathrm{LTE}}/{{\rm{\Lambda }}}_{{{\rm{H}}}_{2},{\rm{n}}\to 0}},\end{eqnarray} \tag{ 37 }$$

where ${{\rm{\Lambda }}}_{{{\rm{H}}}_{2},\mathrm{LTE}}$ is the cooling rate at high densities when ${{\rm{H}}}_{2}$ level populations are in local thermal equilibrium, and ${{\rm{\Lambda }}}_{{{\rm{H}}}_{2},{\rm{n}}\to 0}$ is the cooling rate in the low density limit. We use Equations (6.37) and (6.38) in Hollenbach & McKee (1979) to calculate ${{\rm{\Lambda }}}_{{{\rm{H}}}_{2},\mathrm{LTE}}={{\rm{\Lambda }}}_{{{\rm{H}}}_{2},\mathrm{LTE},\mathrm{rot}}+{{\rm{\Lambda }}}_{{{\rm{H}}}_{2},\mathrm{LTE},\mathrm{vib}}$ for ${{\rm{H}}}_{2}$ rotational and vibrational lines. ${{\rm{\Lambda }}}_{{{\rm{H}}}_{2},{\rm{n}}\to 0}$ is obtained by the fitting functions provided in Glover & Abel (2008), including collisional species H, ${{\rm{H}}}_{2}$, He, ${{\rm{H}}}^{+}$, and e. We use the fitting parameters in their Table 8, assuming a fixed ortho-to-para ratio of ${{\rm{H}}}_{2}$ of 3:1. Glover & Abel (2008) only gives the fit for ${{\rm{H}}}_{2}$ cooling in the temperature range $20\,{\rm{K}}\lt T\lt 6000\,{\rm{K}}$. Similar to the treatment in CO rotational line cooling, we artificially cut off ${{\rm{H}}}_{2}$ cooling below $10\,{\rm{K}}$ and use the ${{\rm{H}}}_{2}$ cooling rate at $T=6000\,{\rm{K}}$ for $T\gt 6000\,{\rm{K}}$. In reality, ${{\rm{H}}}_{2}$ cooling rate is likely to be negligible at $T\lt 10\,{\rm{K}}$.

B.2.4. Dust Thermal Emission

By exchanging energy with dust via collision, gas can either heat or cool. The rate of gas−dust energy exchange is given by (Krumholz 2014)

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{gd}}={\alpha }_{\mathrm{gd}}{{nZ}}_{{\rm{d}}}{T}_{{\rm{g}}}^{1/2}({T}_{{\rm{d}}}-{T}_{{\rm{g}}}),\end{eqnarray} \tag{ 38 }$$

where Z_d is the dust abundance relative to the Milky Way, and ${\alpha }_{\mathrm{gd}}$ the dust−gas coupling coefficient. Positive ${{\rm{\Psi }}}_{\mathrm{gd}}$ means heating of the gas, and negative means cooling of the gas. In the realistic ISM, dust almost always acts to cool the ISM with ${T}_{{\rm{d}}}\lt {T}_{{\rm{g}}}$. We use ${\alpha }_{\mathrm{gd}}=3.2\times {10}^{-34}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{3}\,{{\rm{K}}}^{-3/2}$ for ${{\rm{H}}}_{2}$-dominated regions (Goldsmith 2001). ${\alpha }_{\mathrm{gd}}$ is expected to be a factor ∼3 higher in H-dominated regions (Krumholz et al. 2011), but dust cooling is likely to be only important in dense regions dominated by ${{\rm{H}}}_{2}$.

The dust temperature ${T}_{{\rm{d}}}$ necessary for evaluating ${{\rm{\Psi }}}_{\mathrm{gd}}$ is determined by the total rate of change of dust specific energy per H nucleus (Krumholz 2014):

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}_{{\rm{d}},\mathrm{sp}}}{{dt}}={{\rm{\Gamma }}}_{\mathrm{ISRF}}+{{\rm{\Gamma }}}_{{\rm{d}},\mathrm{line}}+{{\rm{\Gamma }}}_{{\rm{d}},\mathrm{CMB}}+{{\rm{\Gamma }}}_{{\rm{d}},\mathrm{IR}}-{{\rm{\Lambda }}}_{{\rm{d}}}-{{\rm{\Psi }}}_{\mathrm{gd}}.\end{eqnarray} \tag{ 39 }$$

Because the dust specific heat is very small compared to the gas, we can always assume the dust temperature is in equilibrium, ${{de}}_{{\rm{d}},\mathrm{sp}}/{dt}=0$. Note that for very small dust with size a ≲ 10 Å, they are not in thermal equilibrium and can have temperature spikes (Draine 2011). Here we mainly consider bigger grains. In principle, given all the terms in Equation (39), one can solve for ${T}_{{\rm{d}}}$. However, because Equation (39) is a non-linear function of ${T}_{{\rm{d}}}$, this procedure will involve a costly root-finding algorithm. For the purposes of this paper, we assume a fixed dust temperature ${T}_{{\rm{d}}}=10\,{\rm{K}}$. In the density range $n\lt {10}^{4}\,{\mathrm{cm}}^{-3}$ that we consider, dust cooling is never important and contributes to less than ∼5% of the total cooling rate.

B.2.5. Recombination of e on PAHs

B.2.6. Collisional Dissociation of H₂ and Collisional Ionization of H

Here we show the comparison with the PDR code, our chemical network, and the NL99 network. The chemistry is evolved until it reaches equilibrium, and we fix the temperature at $T=20\,{\rm{K}}$. All independent species are shown between densities $n=50\,{\mathrm{cm}}^{-3}$ and $n=1000\,{\mathrm{cm}}^{-3}$, supplementing Figures 2 and 4. Figures 12 and 14 plot the abundances, and Figures 13 and 15 plot the integrated column densities of different species.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We carry out comparisons similar to Appendix C.1, but here we solve the chemistry and temperature simultaneously, until they both reach equilibrium. The results are shown in Figures 16 and 17. Again, our network agrees well with the PDR code, whereas the NL99 network shows significant disagreement. However, the equilibrium temperature in the NL99 network is very similar to that in our network. For example, at $n=100\,{\mathrm{cm}}^{-3}$ and ${A}_{V}\gt 1$, although there is much less CO in the NL99 network, C provides most of the cooling, resulting in similar equilibrium temperatures. As noted by Glover & Clark (2012a), temperature is insensitive to detailed chemistry, although shielding is important because it sets the photoelectric heating rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 18 and 19 plot comparisons at Z = 0.1 between our network and the PDR code. Overall, there is a good agreement. Comparing to Z = 1 cases, CO forms at a much higher density $n\gtrsim 1000\,{\mathrm{cm}}^{-3}$. We note that at $n=1000\,{\mathrm{cm}}^{-3}$, our network overproduces CO in shielded regions. This is partly due to the difference in grain-assisted recombination rates, which is discussed in Appendix C.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Ions such ${{\rm{H}}}^{+}$, ${{\rm{C}}}^{+}$, ${\mathrm{He}}^{+}$, and ${\mathrm{Si}}^{+}$ can recombine when they collide with small dust grains know as PAHs. The recombination rate on dust grains is often higher than direct recombination with electrons, and therefore can have a major effect on chemistry. We use the results of Weingartner & Draine (2001a) for grain-assisted recombination rates in our network (see Table 2). In Weingartner & Draine (2001a), they assume a certain size distribution of PAHs and calculate the charge distribution on PAHs. The final grain-assisted recombination rates are obtained by averaging over the PAH population. The PDR code uses a slightly different approach in Wolfire et al. (2008). The PAHs are assumed to have a single size, and three charge states, ${\mathrm{PAH}}^{+}$, ${\mathrm{PAH}}^{-}$, and ${\mathrm{PAH}}^{0}$. By comparing with observations, Wolfire et al. (2008) calibrate their grain-assisted recombination rate of ${{\rm{C}}}^{+}$ with the parameter ${\phi }_{\mathrm{PAH}}=0.4$.

Figure 20 shows the comparison of the grain-assisted recombination rate of ${{\rm{C}}}^{+}$ between Wolfire et al. (2008) and Weingartner & Draine (2001a). The PAH populations in Wolfire et al. (2008) are assumed to be in charge equilibrium states with reactions described in Appendix C2 of Wolfire et al. (2003), and the total PAH abundance is ${n}_{\mathrm{PAH}}/n=2\times {10}^{-7}$. The reaction rate depends mainly on the parameter $\psi =1.7\chi \sqrt{T}/{n}_{e}$, and only has a very weak dependence on temperature. As can be seen in Figure 20, the Wolfire et al. (2008) and Weingartner & Draine (2001a) rates can differ by a factor of ∼2. In the main part of this paper, we use 0.6 times the Weingartner & Draine (2001a) rates to match with Wolfire et al. (2008) at $\psi \sim {10}^{4}$. Figures 21 and 22 show comparisons with the original Weingartner & Draine (2001a) rates. At $n=100\,{\mathrm{cm}}^{-3}$, The abundances of CO, ${\mathrm{OH}}_{x}$, ${\mathrm{HCO}}^{+}$, and ${\mathrm{He}}^{+}$ can change by a factor of ∼5 around ${A}_{V}\sim 1$, but the chemical abundances at higher density $n=1000\,{\mathrm{cm}}^{-3}$ are largely unaffected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the grain-assisted recombinations, dust grains can also affect the gas-phase chemistry by grain surface reactions (Hollenbach et al. 2009). This includes the formation of OH and H₂O on the surface of dust grains, and the freeze-out of molecules such as CO and H₂O on dust grains in cold and shielded regions. Figures 23 and 24 show the comparisons of the chemical abundances with and without the grain surface reactions. The formation of OH and H₂O on grains results in higher ${\mathrm{OH}}_{x}$ and CO abundances at $n=1000\,{\mathrm{cm}}^{-3}$ and ${A}_{V}\lt 1$. The freeze-out can greatly reduce the abundances of CO and other species at ${A}_{V}\gtrsim 2$. Since the CO abundance is still low in regions that are affected by the formation of OH and H₂O on grains, and the freeze-out depends on the grain size distribution, which is very uncertain in dense and shielded clouds, we do not include grain surface reaction in our network. Ideally, the effect of freeze-out on carbon bearing species should be calibrated with observations of diffuse and dense clouds as was done with oxygen-bearing species in Hollenbach et al. (2009, 2012) and Sonnentrucker et al. (2015).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 25 shows a comparison between our chemical network and the PyPDR code. We disabled the grain-assisted recombination of ions in our network to make a direct comparison with PyPDR, which does not have these reactions. We also updated the reaction rates in the PyPDR code, to be consistent with the reaction rates in this paper (Tables 1 and 2). Note that PyPDR does not include metals such as Si, which is in our network. Figure 25 shows that our chemical network (without grain-assisted recombination) agrees very well with the results from the PyPDR code.

Zoom In Zoom Out Reset image size