WATER, O₂, AND ICE IN MOLECULAR CLOUDS

The numerical code we have developed to model the chemical and thermal structure of an opaque cloud externally illuminated by FUV flux is based on our previous models of PDRs, which are described in Tielens & Hollenbach (1985), Hollenbach et al. (1991), and Kaufman et al. (1999, 2006). This one-dimensional (1D) code modeled a constant density slab of gas, illuminated from one side by an FUV flux of 1.6 × 10⁻³G₀ erg cm⁻² s⁻¹. The unitless parameter G₀ is defined above in such a way that G₀ ∼ 1 corresponds to the average local interstellar radiation field in the FUV band (Habing 1968). The subscript "0" indicates that the flux is the incident FUV flux, as opposed to the attenuated FUV flux deeper into the slab. The code calculated the steady-state chemical abundances and the gas temperature from thermal balance as a function of depth into the cloud. It incorporated ∼45 chemical species, ∼300 chemical reactions, and a large number of heating mechanisms and cooling processes. The chemical reactions included FUV photoionization and photodissociation; cosmic-ray ionization; neutral–neutral, ion–neutral, and electronic recombination reactions, including reactions with charged dust grains; and the formation of H₂ on grain surfaces. No other grain surface chemistry was included in these prior models. The main route to gas-phase H₂O and O₂ through these gas-phase reactions is schematically shown in Figure 1. It starts with the ionization of H₂ by cosmic- rays or X-rays, which eventually leads to H₃O⁺ that recombines to form gas-phase water. It also recombines to form OH that reacts with O to form gas-phase O₂. The code contained a relatively limited carbon chemistry; we included CH and CH₂, as well as their associated ions, but did not include CH₃, CH₄, C₂H, or longer carbon chains, nor any carbon molecules that include S or N.

Zoom In Zoom Out Reset image size

The code modeled regions that lie at hydrogen column densities N ≲ 2 × 10²² cm⁻² (or A_V ≲ 10) from the surface of a cloud. Therefore, it applied not only to the photodissociated surface region, where gas-phase hydrogen and oxygen are nearly entirely atomic and where gas-phase carbon is mostly C⁺, but also to regions deeper into the molecular cloud where hydrogen is in H₂ and carbon is in CO molecules. Even in these molecular regions, the attenuated FUV field can play a significant role in photodissociating H₂O and O₂, and in heating the gas.

The dust opacity is normalized to give the observed visual extinction per unit H nucleus column density that is observed in the diffuse interstellar medium (ISM; e.g., Savage & Mathis 1979; Mathis 1990). The extinction in the FUV, assumed to go as $e^{-b_\lambda A_V}$, is taken from PDR models described most recently in Kaufman et al. (1999, 2006), and mainly derived from the work of Roberge et al. (1991). Here b_λ is the scale factor that expresses the increase in extinction as one moves from the visual to shorter wavelengths. However, the FUV extinction law in molecular clouds is not well known. We discuss in Section 3.2 how our results would vary if b_λ varies.

In the following subsections (Sections 2.2–2.7), we describe the modifications to the basic PDR code that we have made to properly model the gas-phase H₂O and O₂ abundances from the surfaces of molecular clouds to the deep interiors, where the gas-phase elemental oxygen freezes out as water ice mantles on grain surfaces. The modifications include the freeze-out of gas-phase species onto grains, the formation of OH and H₂O as well as CH, CH₂, CH₃, and CH₄ on grain surfaces, and various desorption processes of species from grain surfaces. Figure 2 shows the grain surface formation process for water ice.

Zoom In Zoom Out Reset image size

One key modification is the application of time-dependent chemistry to the deep, opaque interior. Here, certain chemical timescales are comparable to cloud lifetimes, and steady-state chemical networks, such as our PDR code, do not apply. The time-dependent code also increases the number of reactions and species we follow, and this particularly improves the treatment of the carbon chemistry in the cloud interior. Further discussion is found in Sections 2.3, 2.4, and 3.7.

The basic PDR code implicitly assumed an "MRN" (Mathis et al. 1977) grain size distribution n_gr(a) ∝ a^−3.5, where a is the grain radius, extending from very small (≲10 Å) polycyclic aromatic hydrocarbons (PAHs) to large (∼0.25 μm) dust grains. With this distribution, the largest grains provide the bulk of the grain mass, and the intermediate-sized grains (∼100 Å) provide the bulk of the FUV extinction. The smallest particles and PAHs provide the bulk of the surface area, and dominate the grain photoelectric heating process and the formation of H₂ on grain surfaces. The code implicitly assumes an MRN distribution in the sense that the equations for attenuation of incident photons by dust, for grain photoelectric heating, for H₂ formation, and for gas–grain collisions and the charge transfer that may occur in such collisions all assume such a distribution.

For the process of the freeze-out of oxygen species onto grain surfaces, we also implicitly assume an MRN distribution, in the sense that we adopt a cross-sectional grain area σ_H appropriate for an MRN distribution. As will be shown in Section 3, σ_H, which is one-fourth the surface area available for freeze-out, is an important parameter for determining the abundance of gas-phase H₂O in our models. Below in Section 2.5, we show that σ_H is roughly the total cross-sectional area of grains (per H) that are larger than about 20 Å. Smaller grains experience single-photon heating events, which clear them of ice mantles (Leger et al. 1985). In addition, cosmic-rays may clear small grains of ices (Herbst & Cuppen 2006). In an MRN distribution, the cross-sectional area of grains with radii larger than about 20 Å is σ_H ∼ 2 × 10⁻²¹ cm² per H.

We note, however, that the cross-sectional area is somewhat uncertain. It could be smaller if grains coagulate inside of dense molecular clouds. It could be larger due to ice buildup. We note that once all oxygen-based ices freeze out onto such an MRN distribution, it forms a constant mantle thickness of Δa ∼ 50 Å regardless of grain size. While this is insignificant to the large (∼0.1 μm) grains that contain most of the mass, it represents a very large change in radius and area of the smallest (∼20 Å) grains that provide most of the grain surface area. The increased surface area of grains at the end of this freezeout will proportionally increase subsequent gas–grain interaction rates. Thus, neutralization of ions by grains, formation of H₂ on grain surfaces, the freeze-out of species, and the cooling or heating of the gas by gas–grain collisions may all be enhanced by significant amounts in regions where water ice has totally frozen out on grains. However, because of the possibility of grain coagulation, which reduces grain surface area, and because of the uncertainly in the lower size cutoff of effective grains, we have chosen to fix σ_H and treat it as a (constrained) variable with values ranging from σ_H = 6 × 10⁻²² to 6 ×10⁻²¹ cm² per H.

Adsorption is the process by which a gas-phase atom or molecule hits a grain and sticks to the surface due to van der Waals or stronger surface bonding forces. Typically, at the low gas temperatures we will consider, T_gas ≲ 50 K, the sticking probability is of order unity for most species (Burke & Hollenbach 1983). We assume unit sticking coefficients here for all species heavier than helium. Thus, the timescale to adsorb (freeze) is the timescale for a species to strike a grain surface

$$\begin{equation} t_{f,i} \sim [n_{{\rm gr}} \sigma _{{\rm gr}} v_i]^{-1}, \end{equation} \tag{ 1 }$$

where v_i is the thermal speed of species i and where the grain number density n_gr times the grain cross section σ_gr is the average over the size distribution. As discussed above, in an MRN distribution with a lower cutoff size of a ∼ 20 Å, n_grσ_gr ≃ 2 × 10⁻²¹ n cm⁻¹, where n is the gas-phase hydrogen nucleus density in units of cm⁻³. Therefore, t_f,i ∼ 8 × 10⁴(m_O/m_i)^1/2(10⁴ cm^-3/n)(30 K/T)^1/2 years, where m_O is the mass of atomic O and m_i is the mass of the species. Since molecular cloud lifetimes are thought to be 2 × 10⁶ to 2 × 10⁷ years, and the molecules reside in regions with n ≳ 10³ cm⁻³, this short timescale suggests that molecules should mostly freeze out in molecular clouds unless some process desorbs them from grains. Dust vacuums the condensibles quickly.

Consider the depletion of CO in regions where the grain temperatures are too high to directly freeze out CO. CO freezes out when grain temperatures are T_gr ≲ 20 K (see Section 2.4 below). However, H₂O freezes out when grain temperatures are T_gr ≲ 100 K at the densities in molecular clouds. Suppose that 20 K <T_gr < 100 K, so that gas-phase elemental oxygen that is not in CO freezes out as water ice, but CO does not directly freeze. CO is continuously destroyed by He⁺ ions created by cosmic rays, producing C⁺ and O. The resultant O has two basic paths: it can reform CO, or it can freeze out on grain surfaces as water ice, never to return to the gas phase in the absence of desorption mechanisms. The latter process therefore slowly depletes the CO gas-phase abundance, by lowering the available gas-phase elemental O. We see this effect in the late time evolution of CO in Bergin et al. (2000). One effect, not included in Bergin et al., that slows this process is that for grain temperatures ≳20 K, the adsorbed atomic O may be rapidly thermally desorbed before reacting with H, and therefore the production of H₂O ice by grain surface reactions may be substantially curtailed (this will be discussed further in Sections 2.5 and 3.2). Nevertheless, gas-phase H₂O will form and will freeze out, and therefore given enough time, the elemental O not incorporated in refractory (e.g., silicate) grain material is converted to H₂O ice. If elemental C does not evolve to a similarly tightly bound carbonaceous ice, there is the potential for the gas-phase C/O number ratio to exceed unity. This will have consequences for astrochemistry (e.g., Langer et al. 1984; Bergin et al. 1997; Nilsson et al. 2000), and we show some exemplary results in Section 3.5.

Consider as well the case where T_gr ≲ 20 K deep in the cloud (i.e., at A_V > 5; this low T_gr corresponds to G₀ ≲ 500). In this case, the gas-phase CO freezes out with the H₂O, forming a CO/H₂O ice mix. In order for this CO ice to be converted to H₂O ice, the CO needs to be desorbed by cosmic rays. It then reacts in the gas phase with He⁺ as described above, and eventually, the freed O goes on to form mostly water ice. The timescale for this CO cosmic-ray desorption can be quite long, and may dominate the slow evolution of the chemistry toward steady state.

Unfortunately, this timescale is not well known. First of all, cosmic-ray desorption rates are not well constrained. Second, the timescales depend on the structure of the CO/H₂O ice mix, and this structure is not certain. If the molecules are immobile, the CO and H₂O ice molecules are initally mixed and each monolayer contains this mix. However, there is some evidence that either pure CO ice pockets form, or the CO accretes on top of an already formed layer of water ice (Tielens et al. 1991; Lee et al. 2005; Pontoppidan et al. 2004, 2008). Even if the CO and H₂O are initially mixed and immobile, cosmic rays desorb CO much more efficiently than H₂O due to the much smaller binding energy of adsorbed CO compared with adsorbed H₂O (Leger et al. 1985). This could lead to surface monolayers with very little CO compared to the protected layers underneath the surface. The cosmic-ray desorption rate is directly proportional to the fraction of surface sites occupied by CO. In the examples above, this fraction could range from zero to unity. Thus, the desorption timescale could range from large values (≫10⁷ years) to small values (∼10⁵ years). We compute these timescales in Section 2.4.

Because of the moderately long timescale for the depletion of either gas-phase or ice-phase CO via this freezing out of the elemental O in the form of water ice, we apply time-dependent chemical models (see Section 3.7) to the cloud interiors that are run for the presumed lifetimes (∼10⁷ years) of molecular clouds. Time-dependent models are also warranted for clouds of relatively low density, n ≲ 10³ cm⁻³, where t_f ≳ 2 × 10⁶ years (see Equation (1)).

To desorb a species from a grain surface requires overcoming the binding energy that holds the species to the surface. Table 1 lists the binding energies of the key species we treat. Note that water has a very high binding energy and is therefore more resistant to desorption.

Thermal desorption. The rate R_td,i per atom or molecule of thermal desorption from a surface can be written as

$$\begin{equation} R_{{\rm td},i} \simeq \nu _i e^{-E_i/kT_{{\rm gr}}}, \end{equation} \tag{ 2 }$$

where E_i is the adsorption binding energy of species i and $\nu _i= 1.6\times 10^{11}\sqrt{(E_i/k)/(m_i/m_{\rm H})}$ s⁻¹ is the vibrational frequency of the species in the surface potential well, k is Boltzmann's constant, and m_i and m_H are the mass of species i and hydrogen, respectively. One can find the temperature at which a species freezes by equating the flux F_td,i of desorbing molecules from the ice surface to the flux of adsorbing molecules from the gas

$$\begin{equation} F_{{\rm td},i}\equiv N_{s,i} R_{{\rm td},i}f_{s,i} = 0.25 n_i v_i, \end{equation} \tag{ 3 }$$

where N_s,i is the number of adsorption sites per cm² (N_s,i ∼ 10¹⁵ sites cm⁻²), f_s,i is the fraction of the surface adsorption sites that are occupied by species i, n_i is the gas-phase number density of species i, and v_i is its thermal speed. We have assumed sticking probability of unity. From these equations, one can calculate the freezing temperature T_f,i of a species by solving for T_gr

$$\begin{flequation} \fl T_{f,i} \simeq (E_i/k)\left[ \ln \left({{4 N_{s,i}f_{s,i} \nu _i} \over {n_i v_i}}\right) \right]^{-1}, \end{flequation} \tag{ 4 }$$

$$\begin{eqnarray} T_{f,i} &\simeq& (E_i/k)\left[ 57 + \ln \left[ \left({N_{s,i} \over {10^{15}\ {\rm cm^{-2}}}}\right) \left({\nu _i \over {10^{13}\ {\rm s^{-1}}}}\right)\left({{1\ {\rm cm^{- 3}}}\over n_i}\right)\right.\right.\nonumber\\ &&\times\, \left.\left.\left({{10^4\ {\rm cm\ s^{-1}}}\over v_i}\right)\right]\right]^{-1}. \end{eqnarray} \tag{ 5 }$$

This equation shows that at molecular cloud conditions the freezing temperature T_f,i is about 0.02E_i/k, or about 100 K for H₂O but only ∼20 K for CO (see Table 1) We include thermal desorption rates from our surfaces as given by these equations. We compute the number of layers of an adsorbed species on the grain, and only thermally desorb from the surface layer.

Photodesorption. The flux F_pd,i of adsorbed particles leaving a surface due to photodesorption is given by

$$\begin{equation} F_{{\rm pd},i}= Y_i F_{\rm FUV}f_{s,i}, \end{equation} \tag{ 6 }$$

where Y_i is the photodesorption yield for species i, averaged over the FUV wavelength band, and F_FUV is the incident FUV photon flux on the grain surface. We write the FUV photon flux at a depth A_V into the cloud as (see Tielens & Hollenbach 1985)

$$\begin{equation} F_{\rm FUV}= G_0 F_0\, e^{-1.8 A_V}, \end{equation} \tag{ 7 }$$

where F₀ ≃ 10⁸ photons cm⁻² s⁻¹ is approximately the local interstellar flux of 6–13.6 eV photons and G₀ is the commonly used scaling factor (see Section 2.1).

Reliable photodesorption yields are often hard to find in the literature, but for the case of H₂O a laboratory study has been done for the yield created by Lyman α photons (Westley et al. 1995a, 1995b). The mass loss at large photon doses is proportional to the photon flux, showing that desorption is not due to sublimation caused by absorbed heat. The yield is negligible when IR or optical photons are incident. As an approximation, we will assume the average yield in the FUV (912–2000 Å) is the same as the Lyman α (1216 Å) yield. Most of the photodesorption occurs when the incident photon is absorbed in the first two surface monolayers of the water ice (Andersson et al. 2006). The value of the yield (∼10⁻³ − 8 ×  10⁻³) is perhaps ten times smaller than the probability that the incident photon is absorbed in the first two surface monolayers, so the photodesorption efficiency from these layers is of order 10%. We note that this yield from Westley et al. (1995a, 1995b) is similar to the values that K. I. Öberg (2007, private communication) has found in preliminary experiments with a broadband FUV source.

One peculiar result from the Westley et al. (1995a, 1995b) experiments is that the yield first increases with the dose of photons, until it hits a saturation value for that flux. The saturation dose is of order 3 ×10¹⁸ photons cm⁻². In the fields we are considering, G₀∼ 1–10⁵, this dose is achieved in less than 1000 years. Thus, our grains may well be saturated. The interpretation of the increase of the photodesorption yield with dose is that photolysed radicals like O, H, and OH build up on the surface to a saturation value. Chemical reactions of a newly formed radical with a preexisting radical may lead to the desorption of the product (usually H₂O) and explain the dose dependence of the yield. Adding further weight to this interpretation is that the saturated yields decline with decreasing surface temperature ($Y_{\rm H_2O}$ is 0.008 at T_gr = 100 K, 0.004 at 75 K, and 0.0035 at 35 K). As temperatures decrease, the radicals cannot diffuse across the surface as well to react with other radicals, nor can activation bonds be as easily overcome. Since presumably radicals like O, H, and OH are being created, and it is the process of reforming more stable molecules that may help initiate desorption and may kick off neighboring molecules, other molecules may desorb, such as OH, H, O₂, and H₂. These are indeed observed, but Westley et al. (1995a, 1995b) report that most of the desorbing molecules are water molecules. We note, however, a recent theoretical calculation by Andersson et al. (2006) that suggests that one might expect OH and H to be major desorption products, but confirming that the yield of H₂O photodesorption may still be of order 10⁻³. We also find (see Section 4.1) that adopting a yield for OH + H photodesorption from water ice, which is twice that of H₂O photodesorption, provides better agreement to the OH/H₂O abundance ratio measured in translucent clouds. We define Y_OH,w as the yield of OH photodesorbing from water ice; $Y_{{\rm OH},w}=2Y_{\rm H_2O}$. This distinguishes Y_OH,w from Y_OH, which is the yield of OH photodesorbing from a surface covered with adsorbed OH molecules.

The interpretation that radical formation in the surface layer is important in creating photodesorption leads one to re-examine the issue of whether grains near the surfaces of molecular clouds will indeed be saturated with radicals on their surfaces. While it is true that the photon dose is sufficient, this dose was delivered over a long timescale. Westley et al. (1995b) point out that it is possible in this case that radicals may recombine in the time between photon impacts within an area of molecular dimensions (∼10⁻¹⁵ cm²), which might result in smaller yields than the saturated ones.

In light of the uncertainties in the wavelength dependence of the yield over the FUV range, the photodesorption products, and the level of saturation, we adopt $Y_{\rm H_2O} = 10^{-3}$ and Y_OH,w = 2 × 10⁻³ as our "standard" values, but we will show models in which these yields are varied by a factor of 10 from their standard values, maintaining their ratio of 0.5. We will also use the observations of A_V thresholds for water ice formation to constrain the absolute values of these yields (Section 4.2.1), and the observations of OH and H₂O in diffuse clouds to help determine their ratio (Section 4.1).

Recent laboratory measurements of CO photodesorption imply a yield similar to that for H₂O (Öberg et al. 2007). Therefore, we adopt Y_CO = 10⁻³ as our standard value for CO. Unless T_gr < 20 K, adsorbed CO is anticipated to be rare because it thermally desorbs so rapidly. However, in cases with low G₀, CO may occupy ∼1/3 of the surface sites at the H₂O freeze-out column.

For the other species there is less or no reliable photodesorption data in the literature. Since our grain surface is largely H₂O ice, and the photodesorption occurs via the production of radicals (and possibly, the reformation of H₂O from radicals), we assume that the yield for O is Y_O = 10⁻⁴ (a smaller yield since no reformation is possible). The yields for O, OH, and O₂ are less important, because in our models they rarely occupy many of the surface sites. The first two are rare because they rapidly react with H to form H₂O. Table 1 lists the photodesorption yields that we adopted.

A complication arises if grains move rapidly from surface regions (higher temperature and possibly lower density) where H₂O ice can exist to deeper (lower temperature and higher density) regions where CO ice forms and then back again to the surface region. The fraction of the surface now covered by H₂O ice, $f_{s,{\rm H_2O}}$, is now nearly zero, since CO covers the surface. Therefore, H₂O and OH photodesorption rates are very low. We have ignored this effect, assuming a mix of CO and H₂O ice in each adsorbed monolayer.⁷ In this paper, the region where H₂O and OH photodesorption are important is the intermediate or "plateau" region where gas-phase H₂O and OH abundances peak. Here, the photodesorption and/or thermal desorption timescales of CO are much shorter than the dynamical timescales to bring CO-coated ice grains from the deep interior to the plateau region.

Desorption by chemical reaction. There have been some attempts to model the desorption of molecules due to chemical reactions on grain surfaces. D'Hendecourt et al. (1982) and Tielens & Allamandola (1987) discuss an exploding grain hypothesis in which the UV photons produce radicals in the grains, which are relatively immobile until the grain temperature is increased. At this point, they begin to diffuse and react, and the chemical heat causes a chain reaction exploding the grain and "desorbing" species in large bursts. Our photodesorption is a sort of steady-state version of this. O'Neill et al. (2002) examine the effects of hydrogenation on grain surfaces, and in some of their models they assume that the reaction of a newly adsorbed H atom with radicals on the surface, such as OH, will release enough energy to desorb the resultant hydrogenated molecule. In this case, it is the chemical energy carried in by the H atom, and not the photodissociation of the ice itself, which leads to desorption. In our standard case, we assume that H₂ desorbs on formation because of its low adsorption binding energy and its low mass, but that other forming molecules do not. However, we ran one case (see Section 3.5) in which we include this desorption process for all forming molecules in our grain surface chemistry.

Desorption by cosmic rays and X-rays. Deep in the cloud, cosmic-ray desorption of ices becomes important for maintaining trace amounts of O-based and C-based species in the gas phase. Cosmic-ray desorption rates per surface molecule (i.e., molecule in the top monolayer of ice) for various species are listed in Table 1. As discussed above, the desorption rates depend on the abundance of a particular adsorbed species in the surface monolayer. In our standard steady-state models, we assume that if there are a number of species (e.g., CO, CH₄, H₂O) that are adsorbed, that they are well mixed, and each monolayer of ice contains the same mix. However, as discussed above, for CO the fraction of the surface layer containing CO may range from zero to unity. The timescale for the cosmic-ray desorption of a species i is given as

$$\begin{equation} t_{{\rm cr},i} = {{x(i)}\over { 4 R_{{\rm cr},i} N_{s,i} f_{s,i} \sigma _H}}, \end{equation} \tag{ 8 }$$

where x(i) is the abundance of the adsorbed ice species and R_cr,i is the cosmic-ray desorption rate per surface atom or molecule and is given in Table 1. The cosmic-ray desorption time for H₂O_ice is long (∼10¹⁰ years) and is not significant in our models. However, the timescale for CO desorption under certain circumstances may be comparable to the age of the cloud, and can therefore lead to a source of gas-phase CO, which can then serve as a source of gas-phase H₂O and O₂, as discussed above. Substituting from Table 1 for R_cr,CO, assuming an initially high abundance of CO ice, x(CO_ice) = 10⁻⁴, and taking f_s,CO = 0.3, we obtain t_cr,CO ∼ 2.3 × 10⁶ years. However, both R_cr,CO and f_s,CO are uncertain. In the case of the fractional surface coverage f_s,CO, the cosmic-ray desorption process itself, acting on the top monolayer, could result in a very low value of f_s,CO by essentially purging the surface of CO and leaving a fairly pure monolayer of H₂O_ice (see Section 2.3). On the other hand, there may be evidence for a nonpolar CO ice surface on top of a polar H₂O ice mantle resulting in f_s,CO = 1 (see Section 2.3). Because of the uncertainly in the product f_s,COR_cr,CO, we will examine in Section 3.7 a range of values of this parameter extending from 10⁻¹⁵ to 10⁻¹¹ s⁻¹. The corresponding timescale to remove all the ice is about 3 × 10⁵ to 3 × 10⁹ years, ranging from less than to much greater than cloud ages.

Cosmic rays also lead to the production of FUV photons (Gredel et al. 1989) at a level equivalent to G₀ ∼ 10⁻³ throughout the cloud. We include this flux in determining both the photodissociation of gas-phase molecules and the photodesorption of ices.

Following Leger et al. (1985) we assume that, in general, cosmic-ray desorption dominates over X-ray desorption and we therefore ignore this process. Of course, in regions close to X-ray sources and where the UV is shielded by dust, X-ray desorption will become more significant.

Several timescales are needed to fully understand the freezing process. Assuming that the grain size distribution extends to very small grains (or large molecules) whose size is of order a few Angströms, we first compute timescales relevant for single photons to clear grains of ice due to transient heating by the photon. The peak temperature that a silicate grain of radius a achieves after absorbing a typical ∼10 eV FUV photon is given (Drapatz & Michel 1977; Leger et al. 1985)

$$\begin{equation} T_{\rm max}\sim 100 \left({{20\ {\rm \AA}}\over a}\right)^{1.3} \ {\rm K}. \end{equation} \tag{ 9 }$$

The timescale for this grain to radiate away the energy of this photon is

$$\begin{equation} t_r \sim 100\left({{100\ {\rm K}}\over T_{\rm max}}\right)^5 \left({{20\ {\rm \AA}}\over a} \right)^3\ {\rm s}. \end{equation} \tag{ 10 }$$

On the other hand, the timescale for a single H₂O molecule to thermally evaporate from an ice surface on the grain is given by

$$\begin{equation} t_{{\rm evap},{\rm H_2O}} \sim 4 \times 10^{-13}\, e^{{4800\ {\rm K}}\over T_{\rm max}} \ {\rm s}, \end{equation} \tag{ 11 }$$

where ΔE/k = 4800 K is the binding energy of the H₂O molecule (Fraser et al. 2001) and 4 × 10⁻¹³ s is the period of one vibration of the molecule in the ice lattice (see Equation (2) and the following text). For comparison with the above equations, $t_{{\rm evap},{\rm H_2O}}\simeq 10^{29}$ s at T_max = 50 K, 3 × 10⁸ s at T_max = 100 K, and 30 s at T_max = 150 K. Comparison with t_r indicates that the ice mantle will evaporate before cooling when T_max ∼ 150 K. Comparison with the equation for T_max shows that a≳ 15–20 Å is required for grains to be large enough such that single-photon events do not clear ice mantles. Thermal evaporation of H₂O will also cool the grain. Since the binding energy of an H₂O molecule is ∼0.5 eV, ∼20 H₂O molecules must evaporate to offset the heating due to the absorption of a single UV photon. However, we show below that 10 eV photons are absorbed by 20 Å grains more rapidly than O atoms hit the grains to form ice. Therefore, thermal evaporation cannot significantly cool the transiently heated grains compared with radiative cooling, and an entire ice mantle never builds up on the grain surface. The conclusion holds that a ≲ 20 Å grains lose their ice in the presence of single-photon heating events. In summary, for the purposes of modeling oxygen freeze out and ice formation, we assume that the effective surface area for ice freezeout is given by the area of grains that are larger than ∼20 Å.

Several more timescales reveal the fate of an oxygen atom once it sticks to a grain. The timescale t_gr,O for a given grain of radius a to be hit by a gas-phase oxygen atom is given as

$$\begin{eqnarray} &&\hspace{-1.8pc}t_{{\rm gr},{\rm O}} \sim 2.5\times 10^4 \left({10^{-4}\over x({\rm O})}\right)\left({{10\ {\rm K}}\over T_{\rm gas}} \right)^{1/2}\nonumber\\ &&\hspace{-1.8pc}\times\left({{0.1\,\mu \rm m}\over a}\right)^2\left({{10^5\ {\rm cm}^{- 3}}\over n}\right)\ {\rm s}, \end{eqnarray} \tag{ 12 }$$

where x(O) is the gas-phase abundance of atomic oxygen relative to hydrogen. Similarly, the timescale t_gr,H for the same grain to be hit by a hydrogen atom from the gas is given as

$$\begin{eqnarray} \hspace{-1.8pc}t_{{\rm gr},{\rm H}} &&\sim 6.0\times 10^{-1} x(H)^{-1}\left({{10\ {\rm K}}\over T_{\rm gas}} \right)^{1/2}\nonumber\\ \hspace{-1.8pc}\quad\quad &&\times\left({{0.1\,\mu \rm m}\over a}\right)^2\left({{10^5\ {\rm cm}^{-3}}\over n}\right)\ {\rm s}, \end{eqnarray} \tag{ 13 }$$

where x(H) is the abundance of atomic hydrogen. Note that at the column or A_V into the cloud where water begins to freeze out on grain surfaces, x(H) is likely less than unity because most H is in H₂. However, if x(H) ≳ 0.25x(O), then H atoms stick more frequently than O atoms. Provided photodesorption or thermal desorption does not intervene, and assuming that the H atoms are quite mobile on the surface and find the adsorbed O atom, O atoms will react with H atoms on the surface to form OH and then H₂O.⁸

These timescales need to be compared to the timescale for UV photons to photodesorb the O atom and the timescale for O atoms to thermally desorb from grains. The timescale for a grain of size a ≳100 Å to absorb an FUV photon is given as

$$\begin{equation} t_{\gamma } \sim 30 G_0^{-1} e^{1.8A_V} \left({{0.1\,\mu \rm m}\over a}\right)^2 \ {\rm s}. \end{equation} \tag{ 14 }$$

The timescale to photodesorb an O atom from a surface covered with a covering fraction f_s,O of O atoms is then

$$\begin{eqnarray} \hspace{-1.8pc}t_{\gamma -{\rm des}} &&\sim 3\times 10^5 G_0^{-1} e^{1.8A_V} f_{s,O}^{-1}\nonumber\\ \hspace{-1.8pc}\quad\quad &&\times\left({{0.1\,\mu \rm m}\over a}\right)^2 \left({Y_O \over 10^{-4}}\right)^{-1} \ {\rm s}, \end{eqnarray} \tag{ 15 }$$

where Y_O is the photodesorption yield for atomic O. The timescale for thermal desorption of an O atom on a grain is given as

$$\begin{equation} t_{{\rm evap},O} \sim 9\times 10^{-13}\, e^{{800\ {\rm K}}\over T_{\rm max}} \ {\rm s}. \end{equation} \tag{ 16 }$$

Therefore, the timescale for thermal evaporation of an O atom is 5 × 10²² s for a 10 K grain, 2 × 10⁵ s for a 20 K grain, and 0.3 s for a 30 K grain. From Equations (12)–(16) one concludes that if grains are cooler than about 20–25 K and if x(H)>0.25x(O), an adsorbed O atom will react with H to form OH and then H₂O before it is thermally desorbed. Photodesorption of O never appears to be important.

In this paper, we are mainly interested in the abundance of gas-phase H₂O and O₂ at the point where oxygen species start to freeze out in the form of water ice. This "freeze-out depth" A_V ≡ A_Vf is where the gas-phase oxygen abundance plummets and most oxygen is incorporated as water ice on grain surfaces. This occurs when the rate at which relatively undepleted gas-phase oxygen atoms hit and stick to a grain equals the photodesorption rate of water molecules from grain surfaces completely covered in ice, or t_γ−des(A_Vf) ∼ t_gr,O(A_Vf). Using Equations(12) and (15) and taking $Y\equiv Y_{{\rm OH},w} + Y_{\rm H_2O}$ as the total yield of photodesorbing water ice, we find

$$\begin{eqnarray} &&A_{Vf} \simeq 0.56 \ln \left[ 0.8 G_0 \left({Y \over 10^{-3}}\right)\left({{10\ {\rm K}}\over T_{\rm gas}}\right)^{1/2}\left({{10^5\ {\rm cm}^{-3}}\over n}\right)\right],\nonumber\\ \end{eqnarray} \tag{ 17 }$$

if the term in brackets exceeds unity. If it is less than unity, there is freeze-out at the surface of the cloud. As can be seen in Equation (17), freezeout at the surface only occurs for very low values of G₀/n ≲ 10⁻⁵ cm³.

At A_Vf, the timescale t_γf for a grain of size a to absorb an FUV photon can be written by substituting Equation (17) for A_Vf into Equation (14)

$$\begin{eqnarray} &&\hspace{-1.8pc}t_{\gamma f} \sim 24 \left({{10\ {\rm K}}\over T_{\rm gas}}\right) ^{1/2}\left({{0.1\,\mu \rm m}\over a}\right)^2\nonumber\\ &&\hspace{-1.8pc}\quad\times\left({{10^5\ {\rm cm}^{- 3}}\over n}\right)\left({Y \over 10^{-3}}\right)\ {\rm s}. \end{eqnarray} \tag{ 18 }$$

As discussed above, t_γf ≪ t_gr,O, so that at this critical depth no ice forms on small (≲20 Å) grains because single photons can thermally evaporate adsorbed O, OH, or H₂O before thermal radiative emission cools the grain and before another O atom sticks to the grain. Although ice may not form on a ≲ 20 Å grains, they may still contribute to the formation of gas-phase OH or H₂O. Note that once an O atom sticks to a given grain, the timescale t_gr,H for an H atom to stick (Equation (13)) may be shorter than the timescales t_γf for a photon to be absorbed. Thus, the O may be transformed into H₂O on the grain surface before a photon transiently heats the small grain and thermally evaporates the H₂O to the gas phase. Because of these complexities, the possibility of coagulation, and the uncertain nature of small grains, we have chosen to fix the cross-sectional grain area per H nucleus σ_H as constant with depth A_V in a given model, but to then run models with a range of assumed σ_H.

Figure 2 shows the new grain surface oxygen chemistry that we adopt. We maintain the old method of treating H₂ formation on grain surfaces (Kaufman et al. 1999). The main new grain surface chemistry reactions are the surface reactions H + O → OH and H + OH → H₂O. Miyauchi et al. (2008) discuss and review recent laboratory experiments that suggest rapid formation of water ice on grain surfaces. In addition, we include simple grain surface carbon chemistry. Similar to Figure 2, we allow the sticking of C to grains to form C_ice, followed by reactions of surface H atoms to produce CH_ice, CH_2ice, CH_3ice, and CH_4ice, as well as the desorption of all these species. We do not form CO, CO₂, or methanol on grain surfaces in this paper.

In the model, we compare the timescale for O and OH to desorb with the timescale for an H atom to hit and stick to the grain surface. If the desorption timescale is shorter, then we desorb the radical and the grain surface chemistry does not proceed. However, if the desorption time is longer, then we assume that the reaction proceeds. We generally assume that the newly formed surface molecule does not desorb from the heat of formation, but must wait for another desorption process, or, for a further chemical reaction with H, if possible. However, as in O'Neill et al. (2002), we also run a case where the heat of formation desorbs the newly formed molecule. Implicit in this simple grain chemical code is the assumption that the H atom migrates and finds the O or OH before it desorbs or before another H atom adsorbs to possibly form H₂. We note that, in general, it is rare to find even a single O or OH on the grain surface, because desorption or the reaction with H removes them before another O or OH sticks to the surface.

A final change to our chemical code is the inclusion of enhanced rate coefficients for the reactions between atomic or molecular ions with molecules of large dipole moment, such as OH and H₂O. We have used the UMIST rates (e.g., Woodall et al. 2007) and the rates quoted in Adams et al. (1985) and in Herbst & Leung (1986). These rates are larger than the rates without the dipole interactions by a factor ∼10 at T = 300  K and, given their T^−1/2 dependence, increase further at low temperatures.

Most ortho-H₂O detections in ambient (nonshocked) gas are located in Giant Molecular Clouds (GMCs), and, in particular, in the massive cores of these GMCs that are illuminated by nearby OB stars. Ortho-H₂O has not yet been detected toward quiescent "dark clouds" such as Taurus. Therefore, for our standard case, we choose a gas density n = 10⁴ cm⁻³ and an FUV field strength G₀ = 10². We adopt the same gas-phase abundances as in the PDR models of Kaufman et al. (1999) and Kaufman et al. (2006). Details are given in Tables 1 and 2.

In the standard case, the gas (grain) temperature varies from T = 120 K (T_gr = 31 K) at the surface, to T = 22 K (T_gr = 15 K) at the water peak plateau, to T = 16 K (T_gr = 15 K) for A_V = 10. The abundances of O-bearing species in the standard model are shown in Figure 3, while those of C-bearing species are shown in Figure 4. Here, the basic features of the results of freeze-out and desorption are clearly demonstrated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The photodissociation layer. At the surface, the dominant species are precisely those found in gas-phase PDR models. Atomic oxygen is the most abundant O-bearing species, while the transition from C⁺ to CO occurs once significant shielding occurs at A_V ∼ 1.5. Near the surface, all of the oxygen that is not in CO is atomic, and x(O) ⩾ x(CO). At the low dust temperature at the surface, thermal desorption of H₂O is insignificant and the beginnings of an ice layer (less than a monolayer) are evident despite rapid photodesorption by the FUV. In steady state, photodesorption of this ice as well as gas-phase chemistry maintains a small abundance of gas-phase H₂O.

The photodesorbed layer. By A_V ∼ 2.8, enough water has been deposited on grain surfaces to form a monolayer, i.e., $x({\rm H_2O_{\rm ice})}\approx 4 n_{{\rm gr}}\sigma _{{\rm gr}} N_{s,{\rm H_2O}}/n \!\sim\! 8\times 10^{-6} (\sigma _H/2\times 10^{-21}$ cm²). At this point, the water ice abundance increases rapidly with depth. Here, the photodesorption rate decreases with depth due to attenuation of the FUV flux by grains, but the rate of freeze-out only decreases when the gas-phase O abundance decreases. For the gas-phase O abundance to decrease significantly requires many monolayers of ice. This is seen in Figure 3 as the increase in x(H₂O_ice) by two orders of magnitude at A_V ∼ 3. Note the agreement of this freeze-out depth with A_Vf given in Equation (16). As the gas-phase abundance of O decreases, the growth of water ice slows and the ice abundance saturates (i.e., all O nuclei locked in water ice) near A_V ∼ 3.5. Fed by photodesorption, the gas-phase abundance of H₂O maintains a broad plateau x(H₂O) ∼10⁻⁷ from A_V ∼ A_Vf ∼ 3 to A_V ∼ 6, beyond which (i.e., deeper in the cloud) other gas-phase H₂O destruction mechanisms besides photodissociation become dominant. The abundance of gas-phase H₂O near the peak is determined by the balance of photodesorption from the surface layer of ice on the grains with photodissociation of H₂O by FUV photons. Both of these processes are proportional to the local attenuated FUV field; this leads to a plateau in the H₂O abundance independent of G₀ (see Section 3.3).

The abundances of other O-bearing molecules closely track the abundance of H₂O. OH is mainly a product of photodissociation of gas-phase water, while it is destroyed by photodissociation and by reaction with S⁺, both of which depend on the attenuated FUV field. The O₂ is produced primarily through the reaction O + OH → O₂ + H and is destroyed by photodissociation. Thus, the gas-phase H₂O, OH, and O₂ abundances track one another from the beginning of freeze-out through the H₂O plateau. We discuss this in more detail in Section 3.3 below.

The freezeout layer. Beyond A_V > 6, photodesorption of water ice ceases to produce enough gas-phase H₂O or OH to maintain the H₂O and OH plateaus (and therefore the O₂ plateau). Other destruction mechanisms than photodissociation of H₂O take over, and gas-phase abundances of all O-bearing species drop as H₂O ice incorporates essentially all available elemental O.

The destruction of CO by He⁺, followed by freeze-out of the liberated O atom, as described in Section 2.3, is responsible for the rapid drop in the steady-state gas-phase CO abundance at A_V ∼ 6 (see Figure 4) as all of the oxygen becomes locked in water ice. In our rather limited steady-state carbon chemistry, C atoms end up incorporated into CH_4ice once O is liberated from CO. In our time-dependent runs, which have more extensive chemistry, the standard case at A_V = 10 has ∼0.6 of the elemental carbon in CH_4ice and ∼0.4 in CO_ice after 2 × 10⁷ years.⁹ The ratio of the CH_4ice to the H₂O_ice abundance is ∼0.15–0.25 at this time. Öberg et al. (2008) discuss recent observational results on the ice content of grains in molecular clouds with emphasis on CH_4ice; they conclude that CH_4ice does indeed form by hydrogenation on grain surfaces, as we have assumed in our models. However, they observe abundance ratios of CH_4ice to H₂O_ice that are typically 0.05, nearly three to five times smaller than in our standard model. This likely arises because our time-dependent code starts with carbon in atomic or singly ionized form, so that there is more opportunity for C or C⁺ to transform to CH_4ice than if the carbon started in the form of CO. Aikawa et al. (2005) provide theoretical models that result in the observed ratios of CH_4ice/H₂O_ice by initiating carbon in the form of CO. The purpose of this paper is to follow H₂O and O₂ chemistry, and not carbon chemistry. We have therefore also neglected grain surface reactions that could lead to CO_2ice or methanol ice.¹⁰ Our CO_2ice is formed by first forming gas-phase CO₂ by the reaction of gas-phase CO with OH, and then subsequently freezing the CO₂ onto the grain surfaces. In a subsequent paper we plan to study the possible fractionation signature (¹²CO_2ice/¹³CO_2ice) of our gas-phase route to forming CO_2ice in our subsequent paper, which examines carbon chemistry in more detail. In summary, our carbon chemistry here is simplified, and the abundances of carbon ices in the models are very crude.

Dependence on incident FUV flux, G₀. Figure 5 shows the effect of changing the FUV field strength, G₀, on the depth of the gas-phase water peak and the abundance at the peak, while keeping the density fixed. Here G₀ is varied from 1 to 10³. The main effect of changing G₀ is to move the position of the water peak inward with increasing G₀. Varying G₀ from 1 to 10³ changes the depth of freezeout from A_V ∼ 1 to A_V ∼ 7. Gas-phase H₂O peaks once ice has formed at least a monolayer on the grain surfaces, so the location of the water peak also moves inward (logarithmically) with increasing G₀. The peak gas-phase water abundance (∼10⁻⁷) is insensitive to G₀ and insensitive to depth (i.e., A_V) in the water plateau. We also see that the depth where the gas-phase water starts to decline from its peak (plateau) value is also weakly (logarithmically) dependent on G₀, increasing with increasing G₀. Thus, the column of gas-phase water remains fairly constant and independent of G₀.

Zoom In Zoom Out Reset image size

However, some differences arise in the case with G₀ = 10³. Here, the assumption that every O sticking to a grain forms water ice breaks down. The high FUV field absorbed at the cloud surface leads to a high IR field that heats the grains to T_gr ∼ 25 K even at A_V ∼ 8, so that a significant fraction of O atoms is thermally desorbed from grains faster than they can form water ice, resulting in a higher gas-phase abundance of atomic O. At these large depths, the formation of gas-phase water is not mainly by photodesorption of water ice, but by cosmic-ray initiated ion–molecule reactions (see Figure 1). The high gas-phase elemental abundance of O leads to a high rate of formation of H₂O and O₂, whereas the high A_V leads to low photodestruction rates of these two molecules. In addition, although the FUV field is severely attenuated, the residual FUV is sufficient to prevent the complete freeze-out of gas-phase H₂O. The combined effect is a deeper and higher peak in the gas-phase H₂O and O₂ abundances. At very high A_V (A_V > 8), the freeze-out of water ice finally drains the gas-phase oxygen, and gas-phase H₂O and O₂ decrease with increasing A_V. Our time-dependent results at t = 2 × 10⁷ years find most of the oxygen in H₂O-ice (∼70%) and CO₂-ice (∼30%) at A_V > 8, and the carbon is mostly in C₃H-ice (∼60%) and CO₂-ice (∼40%) since the warm grain temperatures lead to thermal desorption of CO-ice and CH₄-ice.

Dependence on gas density, n. In Figure 6, we vary the gas density from 10³ to 10⁵ cm⁻³, keeping the rest of the parameters fixed at their standard values. Comparing this figure with Figure 5, we see that for G₀ < 10³ the plateau water abundance is independent of gas density and G₀, depending only on the photoelectric yield $Y_{\rm H_2O}$ and the grain cross-sectional area per H nucleus σ_H as shown below. The peak abundance is independent of n because both the formation rate per unit volume (photodesorption) and the destruction rate (photodissociation) are proportional to n. For a fixed value of G₀, lowering the gas density moves the water peak further in; the deposition of water ice depends on the product n_grn_O ∝ n², while the photodesorption rate goes only as the grain density n_gr ∝ n times the local FUV field. Therefore, the water peak moves into the cloud with decreasing density. However, the column or A_V where the water starts to decline from its peak (plateau) value is insensitive to n. Thus, the plateau starts to become narrower and narrower (in ΔA_V) as n decreases. As a result, the column of gas-phase water increases with n somewhat. We explain all these trends in more detail below in Section 3.3, where we present a simple analytical analysis.

Zoom In Zoom Out Reset image size

Dependence of gas temperature on G₀ and n. Figure 7 presents the gas temperature as a function of A_V for a range of our G₀, n parameter space. This figure reveals the temperature at which most of the H₂O and O₂ is radiating, as we have marked the peak in the gas-phase water abundance (which coincides with the peak in the O₂ abundance) as solid portions of these curves. We see that over the parameter space G₀ = 1–10⁵, the emission from H₂O and O₂ arises from regions of T∼ 7–50 K, with T generally rising gradually with G₀. For G₀ < 500, there is a true "water plateau" and the solid line corresponds to the plateau region. T∼ 7–30 K in the plateau, and the gas temperature is often higher than the dust temperature because of grain photoelectric heating of the gas. For higher G₀, there is no true plateau but the gas-phase H₂O and O₂ peak at higher A_V, where the gas and dust temperatures are nearly the same. We will show in Section 3.4 that at high A_V, the dust temperature scales roughly as G^0.24₀. Here, we have made solid lines where the H₂O abundance is more than 0.33 of the peak abundance.

Zoom In Zoom Out Reset image size

Dependence on photodesorption yield, $Y_{\rm H_2O}$. Figure 8 presents the profiles of the gas-phase H₂O and O₂ abundances, and the H₂O ice abundance as the photodesorption yield $Y_{\rm H_2O}$ of water ice is varied, maintaining $Y_{{\rm OH},w}=2Y_{{\rm H}_2O}$, but keeping all other parameters standard. In contrast to the lack of dependence of the plateau value of the gas-phase water abundance x_pl(H₂O) on n and G₀, x_pl(H₂O) scales roughly linearly with $Y_{\rm H_2O}$ for low $Y_{\rm H_2O} < 10^{-3}$. This linear dependence can be easily understood: the formation rate (photodesorption) of gas-phase H₂O scales linearly with $Y_{\rm H_2O}$ whereas the destruction rate (photodissociation) at a given A_V is independent of $Y_{\rm H_2O}$.¹¹ On the other hand, at low values of $Y_{\rm H_2O}< 10^{-3}$, the abundance of O₂ rises with $Y_{\rm H_2O}^2$, or, in other words, with [x(H₂O)]².

Zoom In Zoom Out Reset image size

Dependence on effective grain area, σ_H. In Figure 9, we show the effect of varying σ_H with all other parameters at their standard values. The abundance of gas-phase H₂O in the plateau x_pl(H₂O) increases nearly linearly with σ_H. The depth A_V ∼ A_Vf where the H₂O abundance first reaches its plateau value is insensitive to σ_H, as predicted by Equation (17). Similarly, the depth where H₂O finally begins to drop from its plateau value also is insensitive to σ_H. The abundance of O₂, like the abundance of H₂O, scales with σ_H.

Zoom In Zoom Out Reset image size

Dependence on PAHs. We are most interested in the chemistry at A_V ∼ 3–6 where the H₂O and O₂ abundances peak. Because it is not clear if PAHs exist at such depths (Boulanger et al. 1990; Bernard et al. 1993; Ristorcelli et al. 2003), we have not included PAHs in our standard case but have assumed a grain surface area per H given by σ_H, which is appropriate to an MRN distribution extending down to 20 Å grains when computing the neutralization of the gas by collisions with small dust particles. If PAHs are present, the effect of neutralization is enhanced, which changes the gas-phase chemistry somewhat. We have found that the inclusion of PAHs into our standard case does not affect the results appreciably. There is a somewhat greater effect at lower density, such as n ∼ 100 cm⁻³, where the inclusion of PAHs does increase the total water column to high A_V by a factor of about 10, both by broadening the peak plateau and by increasing the peak abundance (for the case G₀ = 1).

Dependence on b_λ. We have run our standard case but have assumed that the small grains responsible for FUV extinction are depleted, for example, by coagulation. In our test run, we have decreased all of our b_λ by a factor of 2, reduced σ_H by 2, and reduced grain photoelectric heating accordingly. The main effect is to increase the characteristic A_V where the H₂O and O₂ plateau by a factor of 2. However, the abundances of H₂O and O₂ in the plateau decrease because of the decrease in σ_H. Therefore, the columns of gas-phase H₂O and O₂ only increase by a factor of about 1.5 from the standard case. The gas temperature does not appreciably change. In summary, the change in the intensities and average abundances of H₂O and O₂ is not large, and is expected to be less than a factor of 2.

Dependence of H₂O and O₂ columns on G₀ and n. In Figure 10, we show the total columns of H₂O and O₂ for a cloud with high A_V > 10 as a function of the incident FUV flux G₀ and the cloud density n. For low G₀ < 500, the H₂O and O₂ columns are roughly independent of G₀ and n, and are of order 10¹⁵ cm⁻² for our standard values of photodesorption yield. These columns arise in the H₂O and O₂ plateaus, at intermediate cloud columns. In the steady-state models, there is so little gas-phase water in the interior of the cloud that the interior does not contribute significantly to the total column. Therefore, once the total column through the cloud N is large enough to include the plateau (N_H ≳ 10²² cm⁻²), the average abundance will scale as N_H⁻¹. Thus, part of the reason for the low average abundances observed by SWAS and Odin is the dilution caused by the lack of H₂O and O₂ at the surface and deep in the interior of very opaque, high-column clouds. At higher G₀ > 500, we see the effect of the thermal evaporation of the O atoms from the grain surfaces. The effect is most dramatic for the gas-phase O₂ column, which rises to values N(O₂) ∼ 2 × 10¹⁶ cm⁻². This is close to what is needed for a detection by Herschel, and therefore it is in this type of environment that Herschel may potentially detect O₂. We note, however, that this result depends on the binding energy of O atoms to water ice surfaces. The value we have adopted (see Table 1) is standard in the literature but it refers to van der Waal binding to a chemically saturated surface. It is quite possible that the O binding energy is larger, and this would then increase the temperature (and G₀) where this effect would be initiated.

Zoom In Zoom Out Reset image size

The results in Sections 3.1 and 3.2 can be understood by a simple analytical chemical model that incorporates the main physics. Such a model, though approximate, has the advantage of allowing one to determine and understand the sensitivity to various model parameters, and serves to validate the numerical model. In this subsection, the chemical rates are all taken from the UMIST database, except for the ion–dipole rates and the photodesorption rates discussed in Section 2.7. In the following, the reader is referred to Table 3 in the Appendix for quick reference to symbols.

The abundance x(H₂O) of gas-phase H₂O can be roughly traced and understood by simple analytical formulae. Photodissociation dominates the destruction of gas-phase H₂O from the surface to the depth A_Vd where gas-phase destruction with H⁺₃ generally takes over. From 0 < A_V < A_Vd, the simple model assumes that H₂O is formed by photodesorption of H₂O ice. Therefore, by equating formation to destruction

$$\begin{equation} G_0 F_0\, e^{-1.8A_V} Y_{\rm H_2O}f_{s, {\rm H_2O}} n_{{\rm gr}}\sigma _{{\rm gr}} = G_0 R_{\rm H_2O}\,e^{-1.7A_V} n({\rm H_2O}), \end{equation} \tag{ 19 }$$

or

$$\begin{equation} x({\rm H_2O}) = F_0 Y_{\rm H_2O} f_{s, {\rm H_2O}} \sigma _H R_{\rm H_2O}^{-1}\,e^{-0.1A_V}, \end{equation} \tag{ 20 }$$

where $R_{\rm H_2O} = 5.9 \times 10^{-10}$ s⁻¹ is the unshielded PDR of H₂O in an FUV field of G₀ = 1. The difference in the b_λ factors of 1.7 and 1.8 is because the two processes are dominated by somewhat different FUV wavelengths. The fractional coverage of the surface by water ice $f_{s,{\rm H_2O}}$ is given by equating the sticking of O atoms to grain surfaces to the photodesorption rate of H₂O (including both desorption of H₂O and of OH + H). Therefore, recalling that $Y= Y_{\rm H_2O} + Y_{{\rm OH},w}$,

$$\begin{equation} n(O)v_On_{{\rm gr}}\sigma _{{\rm gr}} = G_0F_0\,e^{-1.8A_V}Yf_{s,{\rm H_2O}}n_{{\rm gr}}\sigma _{{\rm gr}} \end{equation} \tag{ 21 }$$

or

$$\begin{equation} f_{s,{\rm H_2O}} = {{n(O)v_O\,e^{1.8A_V}}\over {G_0 F_0 Y}}. \end{equation} \tag{ 22 }$$

Once a monolayer forms, $f_{s,{\rm H_2O}} \sim 1$. We have found in our numerical runs that sometimes $f_{s,{\rm H_2O}}$ saturates at ∼0.5 rather than 1.0, because of the presence of other ices such as CO ice or CH₄ ice. However, assuming a saturation value of unity, the critical depth A_Vf at which a monolayer forms is then given as (this is equivalent to Equation (17) in Section 2.5)

$$\begin{equation} A_{Vf}= 0.56 \ln \left({{G_0 F_0 Y}\over {n(O)v_O}}\right), \end{equation} \tag{ 23 }$$

or a critical H nucleus column (Spitzer 1978) from the surface of

$$\begin{equation} N_f = 2\times 10^{21} A_{Vf} \ {\rm cm}^{-2}. \end{equation} \tag{ 24 }$$

We therefore find, using Equations (20) and (22)

$$\begin{equation} x({\rm H_2O}) = {{n(O) v_O Y_{\rm H_2O}\sigma _H\, e^{1.7A_V}} \over {YG_0 R_{\rm H_2O}}} \quad {\rm for} \quad A_V < A_{Vf}, \end{equation} \tag{ 25 }$$

$$\begin{eqnarray} x_{\rm pl}({\rm H_2O}) &=& {{F_0 Y_{\rm H_2O} \sigma _H e^{-0.1A_V}} \over {R_{\rm H_2O}}}\nonumber\\ &\simeq& 2 \times 10^{-7} \left({Y_{\rm H_2O} \over {10^{-3}}}\right)\left({\sigma _H \over {2\times 10^{-21} \ \rm cm^2}}\right)e^{-0.1A_V}\nonumber\\ &&\ \ \quad {\rm for}\quad A_{Vf} < A_V < A_{Vd}. \end{eqnarray} \tag{ 26 }$$

The abundance of gas-phase H₂O rises exponentially with depth (as $e^{1.7A_V}$) until a monolayer of water ice has formed on the grains. It then forms a plateau with value x_pl(H₂O) that is independent of G₀ and n and weakly dependent on A_V. The plateau abundance goes linearly with Y and σ_H. The lack of dependence on G₀ can be understood since both the formation and destruction of gas-phase H₂O depend linearly on G₀. Higher fluxes do decrease the surface abundance of gas-phase water, but the freeze-out depth is deeper, so there is more dust attenuation, and the gas-phase water abundance rises and peaks at the same constant value. More precisely, regardless of the incident FUV flux, the local attenuated FUV flux at the freeze-out point is always the same, because it must be the flux that photodesorbs an ice-covered grain at the rate that oxygen from the gas resupplies the surface with ice. Therefore, the gas-phase water abundance is the same as well, regardless of G₀, only depending on the assumed values of the grain parameters and the photodesorption yield. The lack of dependence on n can be understood similarly since both the formation rate per unit volume and the destruction rate per unit volume depend linearly on n. On the other hand, x(H₂O) ∝ Yσ_H, since the formation rate by photodesorption goes as this product whereas the destruction rate is independent of them.

In order to determine A_Vd and to provide an approximation for x(H₂O) that applies for A_V > A_Vd, we need to include other destruction routes for gas-phase H₂O that begin to dominate at high A_V. We find that generally the route that takes over from photodissociation is reaction with H⁺₃ to form H₃O⁺. H₃O⁺ recombines one-fourth of the time to reform water (not a net destruction) and three-fourths of the time to form OH (Jensen et al. 2000; Neufeld et al. 2002). The net rate coefficient for the destruction of gas-phase H₂O by H⁺₃ is then

$$\begin{equation} \gamma _{H_3^+} = 0.75\big(4.5 \times 10^{-9} T_{300}^{-1/2}\big) \ {\rm cm^3\ s^{-1}} \end{equation} \tag{ 27 }$$

where T₃₀₀ = T/300 K. Equating the formation of gas-phase H₂O by photodesorption to the destruction of H₂O by both photodissociation and H⁺₃, we obtain a more general solution for x(H₂O) for A_V > A_Vf that applies beyond A_Vd

$$\begin{equation} x({\rm H_2O}) \simeq {{x_{\rm pl}({\rm H_2O})} \over {1 + \big[{{\gamma _{H_3^+} n({\rm H_3^+})}\over {G_0 R_{\rm H_2O}}}\big]e^{1.7A_V}}}, \ {\rm} \ \quad A_V > A_{Vf}. \end{equation} \tag{ 28 }$$

Note that the second term in the denominator, which represents destruction by H⁺₃, begins to dominate when its value reaches and exceeds unity. At this point, x(H₂O) begins to drop from its plateau value.

H⁺₃ is formed by cosmic rays ionizing H₂ followed by rapid reaction of H⁺₂ with H₂. When gas-phase CO is depleted, it is destroyed by dissociative recombination with electrons. The electrons are provided by Si⁺ and S⁺ at moderate A_V and by H⁺₃, He⁺, and H⁺ at high A_V, so the electron density is complicated. At high A_V, the code results suggest that the electron density is roughly n(e) ≃ 2n(H⁺₃). We then obtain

$$\begin{equation} n({\rm H_3^+})_{{\rm high}A_V} \simeq 10^{-5}T_{300}^{1/4} n^{1/2}\ {\rm cm^{-3}}, \end{equation} \tag{ 29 }$$

with H nucleus density n in cm⁻³. In Equation (29), we have assumed a total (including secondary ionizations) cosmic-ray ionization rate of 5 × 10⁻¹⁷ s⁻¹ for H₂ (Dalgarno 2006). At moderate A_V, Si⁺ usually provides most of the electrons. Gas-phase Si⁺ is produced by photoionization of gas-phase Si and destroyed by recombination with electrons. The Si⁺ density, and therefore the electrons provided by Si⁺, exponentially drops with A_V as the FUV photons that photoionize Si drop with increasing dust attenuation. As long as most of the gas-phase silicon is Si⁺, the electron density is high and this suppresses the abundance of H⁺₃, prolonging the water plateau to higher A_V (see Equation (28)). However, we find that once the gas-phase Si⁺ density drops to $2n({\rm H_3^+})_{{\rm high}A_V}$, so that the electrons provided by Si⁺ are comparable to those provided by H⁺₃, He⁺, and H⁺, n(H⁺₃) has risen to ${\sim} n({\rm H_3^+})_{{\rm high}A_V}$ and the H⁺₃ term in the denominator of Equation (28) dominates. Here, the gas-phase water abundance starts to fall. This occurs when

$$\begin{equation} A_{Vd} = 0.59 \ln \left[ 4.5 \times 10^5 G_0 n^{-1/2} T_{300}^{-0.2} \left({n_{{\rm Si}T} \over n_{{\rm H_2O}-{\rm gr}}}\right)^{1/2} \right], \end{equation} \tag{ 30 }$$

where n_SiT = 1.7 × 10⁻⁶n is the density of silicon both in the gas phase and as ice mantles on grains, and $n_{\rm H_2O-gr} \sim 3 \times 10^{-4}n$ is the density of water molecules incorporated in ice mantles.

We have checked our analytical expressions for A_Vf (Equation (23)), A_Vd (Equation (30)), and x_pl(H₂O) (Equation (26)) and find good agreement (A_Vf and A_Vd better than ±25% and x_pl(H₂O) good to a factor of 2 over our (G₀, n) parameter space for G₀ ≲ 500). For G₀ > 500, the grains are so warm that O atoms adsorbed to grains thermally evaporate before reacting with H atoms to form OH on the grain surface. As discussed in Section 3.2, this results in a change in the behavior of gas-phase H₂O.

The abundance of O, OH, and O₂ in the H₂O plateau region A_Vf < A_V < A_Vd can also be analytically estimated from the above expression for x_pl(H₂O). We first find $x_{\rm pl}(\it OH)$ by equating formation by photodissociation of H₂O and photodesorption of OH from water ice with destruction by photodissociation of OH and reaction with S⁺ and He⁺. The photodesorption rate of OH is two times the rate of H₂O photodesorption, and therefore two times the rate of H₂O photodissociation. On the other hand, the destruction of OH by S⁺ and He⁺ is roughly two times that of OH photodissociation. Therefore,

$$\begin{eqnarray} &&\hspace{-1.8pc}x_{\rm pl}(OH) \sim \left({{R_{{\rm H_2O}}} \over {R_{\rm OH}}}\right) x_{\rm pl}({\rm H_2O}) \sim 1.7x_{\rm pl}({\rm H_2O}),\nonumber\\ &&\hspace{-1.8pc}\quad A_{Vf} < A_V < A_{Vd} \end{eqnarray} \tag{ 31 }$$

where R_OH = 3.5 × 10⁻¹⁰ s⁻¹ is the unshielded PDR of OH when G₀ = 1. Thus, OH tracks H₂O with somewhat greater abundance in the plateau, as observed in Figure 3.

The gas-phase O abundance is determined by balancing the formation of O by photodissociation of OH with the removal of O by adsorption onto grains

$$\begin{eqnarray} x_{\rm pl}(O) \!\!&\simeq&\!\! \left[{{G_0 R_{\rm OH} e^{-1.7 A_V}}\over {\sigma _{\rm H} n v_O}}\right] x_{\rm pl}(OH)\nonumber\\ \!\!&\simeq&\!\! 4\times 10^{-5} \left({{G_0/n}\over 10^{-2}}\right) \left({e^{-1.7A_V}\over 10^{-3}}\right) \left({x_{\rm pl}({\rm H_2O}) \over 10^{-7}}\right). \end{eqnarray} \tag{ 32 }$$

Therefore, x_pl(O) exponentially declines with $e^{-1.7A_V}$ due to the declining photodissociation of OH, but with a constant (with A_V) collision rate of O with grains, as seen in Figure 3. It also increases with G₀/n for fixed A_V. This formula breaks down when the right-hand side exceeds ∼10⁻⁴, where x_pl(O) saturates since all elemental gas-phase O is now in atomic form.

The O₂ abundance follows from the abundances of O and OH since O₂ is formed by the neutral reaction of OH with O (rate coefficient $\gamma _{O_2}= 4.3 \times 10^{-11}(300\ {\rm K}/T)^{1/2} e^{-30\ {\rm K}/T}$ cm³ s⁻¹). The O₂, like H₂O and OH, is destroyed by photodissociation in the plateau region. Therefore,

$$\begin{eqnarray} x_{\rm pl}(O_2) &\simeq & \left[{{1.7^2 \gamma _{\rm O_2}}\over {\sigma _{\rm H} v_O}}\right] \left({ R_{\rm OH} \over R_{\rm O_2}}\right) e^{0.1A_V} [x_{\rm pl}({\rm H_2O})]^2\nonumber\\ &\sim & 1 \times 10^{-7} e^{0.1A_V}\left({{2\times 10^{-21}\ {\rm cm^2}} \over {\sigma _{\rm H}}}\right)\nonumber\\ &\times& \left({30 \over {T}}\right)e^{-30/T} \left[{x_{\rm pl}({\rm H_2O}) \over 10^{-7}}\right]^2, \end{eqnarray} \tag{ 33 }$$

where $R_{\rm O_2}=6.9\times 10^{-10}$ s⁻¹ is the unshielded rate of photodissociation of O₂ in a G₀ = 1 field. The [x_pl(H₂O)]² dependence, which translates to a $Y_{{\rm H_2O}}^2$ dependence in the plateau for low $Y_{{\rm H_2O}} < 10^{-3}$, is clearly seen in Figure 8. Note also that Equation (33) explains the lack of dependence of x_pl(O₂) on n and G₀, and also predicts that x_pl(O₂) ∝ σ_H as observed in the numerical models.

In summary of the above, the approximate analytical equations can be used in the cases with G₀ ≲ 500 to understand and quantitatively estimate: (1) x_pl(H₂O, OH, O and O₂) and their dependence (or independence) on n, G₀, Y, and σ_H; (2) the depth A_Vf of the onset of the plateau, the depth A_Vd of the termination of the plateau, the width ΔA_V = A_Vd − A_Vf of the plateau, and their dependence on the above parameters. At high $Y_{{\rm H_2O}}> 10^{-3}$, the abundances of O, OH, O₂, and H₂O are overestimated by factors of order 2 in the above formulae because we have ignored the presence of other ices, which reduce the fractional surface coverage of water ice, and therefore the production of gas-phase H₂O and its photodissociation products.

Finally, the above analytical expressions can be used to predict the optimum conditions for observing H₂O and O₂. In general, for our model runs, the intensity of the H₂O ground state ortho line at 557 GHz and para line at 1113 GHz can be estimated by assuming that the emission is "effectively thin."¹² In other words, although the optical depth in the ground state may be larger than unity, the gas density is so far below the critical density (∼10⁸ cm⁻³) that self-absorbed photons are re-emitted and eventually escape the cloud (see discussion in Section 5). In this limit, every collisional excitation results in an escaping photon, and

$$\begin{equation} I_{557} \simeq {{ h\nu [x_o n \gamma _{{\rm ex},o}+x_p n \gamma _{{\rm ex},p}] N(o-{\rm H_2O})}\over {4 \pi }}, \end{equation} \tag{ 34 }$$

where γ_ex,o ≃ 2.5 × 10⁻¹⁰e^−26.7/T and γ_ex,p ≃ 3.5 × 10⁻¹¹e^−26.7/T cm³ s⁻¹ (valid for T ≲ 30 K) are the rate coefficients for the collisional excitation of the 557 GHz transition by impact with ortho- and para-H₂, respectively (Faure et al. 2007; Dubernet et al. 2006). Here, x_o(x_p) is the abundance of ortho(para)-H₂ with respect to H nuclei,¹³ hν is the energy of the 557 GHz photon, and N(o − H₂O) is the column of ortho water. Similarly, the para H₂O 1113 GHz line intensity is

$$\begin{equation} I_{1113} \simeq {{ h\nu [x_o n \gamma _{{\rm ex},o}+x_p n \gamma _{{\rm ex},p}] N(p-{\rm H_2O})}\over {4 \pi }}, \end{equation} \tag{ 35 }$$

where γ_ex,o ≃ 3 × 10⁻¹⁰e^−53.4/T and γ_ex,p ≃ 9.0 ×  10⁻¹¹e^−53.4/T cm³ s⁻¹ (valid for T ≲ 30 K) are the rate coefficients for the collisional excitation of the 1113 GHz transition by impact with ortho- and para-H₂. From the above equations, assuming most of the H₂O column arises in the plateau between A_Vf and A_Vd, the ortho(para)-H₂O column is given as

$$\begin{eqnarray} N(o(p)-{\rm H_2O}) &\simeq& 2\times\ 10^{21}(A_{Vd}-A_{Vf})\nonumber\\ &&\times\ x_{\rm pl}({\rm H_2O}) f_o(f_p)\ {\rm cm^{-2}}, \end{eqnarray} \tag{ 36 }$$

where f_o(f_p) is the fraction of H₂O in the ortho(para) state. Substituting our analytic equations for these parameters, we obtain

$$\begin{eqnarray} &&\hspace{-1.8pc}I_{557} \simeq 3 \times 10^{-7}f_o[x_o +0.14 x_p]\left({n \over { 10^4 \ \rm cm^{-3}}}\right)\nonumber\\ &&\times\left({Y_{{\rm H_2O}} \over {10^{-3}}}\right) \hspace{-1.8pc}\qquad\left({\sigma _H \over {2 \times 10^{-21} \rm cm^2}}\right)\nonumber\\ &&\times\Delta A_V e^{-0.1A_V} e^{-26.7/T} \ {\rm erg\ cm^{-2}\ s^{-1}\ sr^{-1}} \end{eqnarray} \tag{ 37 }$$

and

$$\begin{equation} I_{1113} \simeq 2.4 \left({f_p\over f_o}\right)\left({{x_o + 0.3x_p}\over {x_o + 0.14x_p}}\right) e^{-26.7/T} I_{557}, \end{equation} \tag{ 38 }$$

where ΔA_V = A_Vd − A_vf. Note that the column of H₂O is relatively independent of G₀ and n. However, the intensity I₅₅₇ or I₁₁₁₃ is proportional to the H₂O column times γ_exn ∝ ne^−26.7/T or ne^−53.4/T, and therefore is somewhat sensitive to G₀ (which determines T) but is more sensitive to n. The 557 GHz line depends on f_o, which is a sensitive function of T_gr for low G₀ (i.e., low T_gr), as will be described below. Therefore, the 557 or 1113 GHz H₂O line will be strongest in dense regions illuminated by strong FUV fluxes (assuming sufficient column or A_V to encompass the plateau, and assuming the regions fill the beam). We note that the gas and grain temperatures vary somewhat through the water plateau region. We take the gas T and the grain T_gr to be the average of the values at the freeze-out point A_Vf and at the A_V where the gas-phase water abundance peaks. Figure 11 shows a contour plot of this average gas temperature T_ave as a function of G₀ and n. Over nearly the entire range of G₀ and n, T_ave only varies from ∼ 15–40 K. Recall that Equations (37) and (38) break down for G₀ ≳ 500, but note that from our numerical results I₅₅₇ and I₁₁₁₃ will increase further, and the analytical equations will underestimate the intensities with increasing G₀ above G₀ ∼ 500, since the water column rises here.

Zoom In Zoom Out Reset image size

Similarly, we can estimate the intensity of the O₂ transition. Here, the density is above the critical density (∼10³ cm⁻³) for the 487 GHz O₂ transition, so the levels are in LTE. In addition, because the transition has a low Einstein A value, large columns are needed to make the transition optically thick. Thus, again assuming optically thin emission, the intensity I₄₈₇ ∝ N(O₂)e^−25K/T/Z(T), where Z(T) is the partition function. An approximation to Z(T) for 0 K <T < 25 K is given as Z(T) ≃ 3 + 2.1T + 0.1T². Using the equations from this subsection, we then estimate

$$\begin{eqnarray} &&\hspace{-1.8pc}I_{487} \simeq 6 \times 10^{-8}\left({{\Delta A_V} \over {Z(T)}}\right) \left({{30\ \rm K}\over T}\right) \left({Y_{{\rm H_2O}}\over {10^{-3}}}\right)^2\nonumber\\ &&\hspace{-1.8pc}\quad\times\, \left({\sigma _H \over {2\times 10^{-21} \rm cm^2}}\right) e^{-55/T} e^{-0.05(A_{Vf}+A_{Vd})}\, {\rm erg\ cm^{-2}\ s^{-1}\ sr^{-1}}.\nonumber\\ \end{eqnarray} \tag{ 39 }$$

Note that again, this equation breaks down for G₀ ≳ 500, and the equation will underestimate the intensity for higher G₀ since the O₂ column rises. In general, however, the equation predicts little n dependence (there is a small n dependence in the A_V terms). At low G₀ ≲ 500, the equation predicts only moderate G₀ dependence. Although various terms depend on T, the temperature of the O₂ plateau region is not sensitive to G₀ since the plateau retreats further into the shielded region for higher G₀. In addition, the T dependence in the exp(-55/T) term is partially cancelled by the TZ(T) term in the denominator.

Using the numerical models, Figures 12 and 13 present H₂O I₅₅₇ and I₁₁₁₃ and O₂ I₄₈₇ contours as functions of n and G₀ for fixed standard values of the other parameters. We have assumed that the H₂O ortho-to-para ratios are in LTE with the grain temperature (which in the plateau is nearly the same as the gas temperature except for very low ratios of G₀/n, which bring the plateau very close to the surface of the cloud where T > T_gr). The grain temperature is appropriate since a given H₂O molecule forms on a grain and spends too little time in the gas to equilibrate to the gas temperature before it is photodissociated. The grain temperature in the plateau is approximately (Hollenbach et al. 1991)

$$\begin{equation} T_{{\rm gr},{\rm pl}} \simeq \big(3\times 10^4 + 2 \times 10^3 G_0^{1.2}\big)^{0.2}. \end{equation} \tag{ 40 }$$

A fit to the LTE ortho-to-para ratio of H₂O valid to 10% in the region 5 K <T_gr < 20 K and good to 20% for T_gr > 20 K (e.g., Mumma et al. 1987) is given as

$$\begin{equation} {f_o\over f_p} = 0.8 e^{-{{22\,{\rm K}} \over T_{{\rm gr}}}} + 2.2\, e^{-\left({{15\,{\rm K}} \over T_{{\rm gr}}}\right)^2}. \end{equation} \tag{ 41 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The H₂ ortho-to-para ratio is calculated in steady state using the rates and processes described by Burton et al. (1992). They are neither 3 nor LTE. They are insensitive to n in the plateau, but sensitive to G₀, which affects T. The H₂ ortho-to-para ratio x_o/x_p ≃ 0.5 at G₀ = 3 × 10⁴, 0.3 at G₀ = 1.5 × 10⁴, and 0.1 for G₀ = 4 × 10³. For lower values of G₀, x_o/x_p < 0.1 so that para H₂ collisions dominate the excitation of water.

The 557 and 1113 GHz intensities behave as described by the analytical formulae (Equations (37) and (38)) presented in Section 3.3 for G₀ ≲ 500. The H₂O intensities rise roughly linearly with n. For a fixed density, I₅₅₇ and I₁₁₁₃ tend to rise with G₀ (which raises T_gr as shown in Equation(40), and which tends to raise T as can be seen in Figure 11) partly due to the excitation factors e^−26.7/T and e^−53.4/T, but also because increasing T increases the ortho/para ratio of H₂, and the ortho excitation rates in collisions with H₂O are much higher. I₁₁₁₃ does not rise as much with increasing G₀ as one might naively expect given its excitation factor e^−53.4/T, because the increasing T reduces f_p from a value of ∼1 at low G₀ to 0.25 at high G₀.¹⁴

The following summarizes the accuracy of the analytical approximations for the H₂O and O₂ intensities for $Y_{{\rm H_2O}} \sim 10^{-3}$. For G₀ < 500, the analytical expression for I₅₅₇ agrees with the numerical results to within a factor of 3 for n = 10³–10⁵ cm⁻³. The agreement for I₁₁₁₃ is within a factor of 3 for n = 10⁴–10⁵ cm⁻³, but only to a factor of 5 (with the analytical underpredicting the intensity) for n = 10³ cm⁻³. However, for this low value of n, the intensities are weak and not detectable. For higher values of G₀, the analytical formula underpredicts the intensity because of the larger columns of H₂O and O₂ (see Figure 10 and Section 5). For higher values of n, the water plateau is so close to the surface that the gas temperature drops rapidly and significantly from A_Vf to A_Vd; here, the analytical formula is not accurate since it assumes constant average temperatures for both grains and gas. Nevertheless, the analytical formulae are useful in making estimates of line intensities, and also illustrate the dependences of the intensities on the physical parameters.

Following O'Neill et al. (2002), as discussed in Section 2.4, we have run our standard case with the modification that all molecules formed on grain surfaces instantly desorb. In our standard case, only H₂ instantly desorbs. The basic differences are that the freeze-out depth A_Vf moves inward (increases) by about two magnitudes of visual extinction due to the chemical release of adsorbed OH and H₂O, the peak gas-phase abundance of H₂O increases by about a factor of 10, and that of O₂ by a factor of roughly 100. We do not pursue this type of model further here, but note that it could warrant further study, although it appears to predict excessive H₂O and O₂ intensities. In the original O'Neill et al. work, photodesorption and cosmic-ray desorption were not included.

In order to test the importance of the formation of OH, H₂O, and carbon molecules on grain surfaces in our model, we have run the standard case (G₀ = 100, n = 10⁴) and two other cases (G₀ = 100, n = 10³ cm⁻³ and G₀ = 10, n = 10⁴ cm⁻³) with only H₂ formed on grain surfaces and no other grain chemistry. Like the model described above in Section 3.5, one main effect is to push the peak in the H₂O and O₂ abundances to higher A_V, since grain formation of H₂O is suppressed. However, in this case, there is no enhancement of the peak abundances because the formation rates of OH and H₂O are suppressed even more than in the case studied in Section 3.5. In addition, there can be significant depletion of gas-phase elemental O at low G₀ by the formation of atomic O ice, which does not occur in Section 3.5 or in our models with grain chemistry in Sections 3.1 and 3.2. In the (G₀ = 100, n = 10⁴ cm⁻³) case and the (G₀ = 10, n = 10⁴ cm⁻³) case, there is a reduction in both the columns of H₂O and O₂ and in the gas and grain temperatures of the emitting region compared to the cases with grain chemistry. Thus, the predicted lines are significantly weaker. Because of the good quality of our grain chemistry model fits to the data described in Section 4, and because of the laboratory and observational evidence for grain chemistry described in Sections 2.3 and 3.1, we do not consider this case further.

The numerical results and analytical solutions presented so far are for steady-state solutions to the chemistry. However, as noted in Sections 2.3 and 2.4, in the opaque interior of molecular clouds at A_V > A_Vd, the chemical timescales may exceed or be comparable to the cloud age, so that time-dependent solutions to the chemistry are needed. The long chemical timescales include the freeze-out of gas species onto grains for n ≲ 10³ cm⁻³ (see Equation (1)), the cosmic-ray desorption of CO ice (see Section 2.4), and the conversion of gas-phase CO to water ice, which is initiated by He⁺ (see Section 2.3). In the opaque "freeze-out" regions, the attenuated FUV flux plays little role in the chemistry so that the time-dependent abundances do not vary with A_V. Therefore, we present results at a representative A_V = 10 for the abundances of gas and ice species as a function of the age of the cloud, the density of the cloud, and the cosmic-ray desorption rate of CO.

The time-dependent code is that of Bergin et al. (2000), updated to have the same grain chemistry and desorption processes as described above for the steady-state code. It includes more species than the PDR surface code, and therefore tracks the carbon chemistry more accurately in the inner regions. It takes as initial conditions the chemistry of translucent clouds: all of the hydrogen is initially molecular H₂ with the other species atomic (for ionization potentials >13.6 eV) or singly ionized otherwise.

The results are shown in Figures 14 and 15, where we present the abundances of O₂, H₂O, and CO plotted as a function of the cosmic-ray desorption rate of CO ice, for two times representative of cloud ages (2 × 10⁶ years in Figure 14 and 2 × 10⁷ years in Figure 15), for the standard G₀ = 100, and for three densities (10³, 10⁴, and 10⁵ cm⁻³). Of particular note is Figure 15, where one sees the gas-phase molecules increase in abundance as the cosmic-ray desorption rate increases, reach a peak when the timescale for cosmic-ray desorption of CO is about 0.1 − 1.0 times the cloud evolution time of 2 × 10⁷ years, and then decrease as the cosmic-ray rate increases further and the desorption time becomes less than ∼2 × 10⁶ years. If the rate is too low, then the CO ice remains on the grain and the gas-phases abundances are low. If the rate is too high, then the CO ice rapidly desorbs, but then the CO in the gas breaks down by reaction with He⁺ and the gas-phase elemental oxygen freezes out as water ice, thereby lowering the gas-phase abundances of CO, H₂O, and O₂. The main point is that in all cases for t = 2 × 10⁷ years and for most cases with t = 2 × 10⁶ years, the abundances of H₂O and O₂ at high A_V, even when time-dependent chemistry is taken into account (steady-state chemistry is more extreme in this respect as seen in Figure 3), are substantially lower than the peak abundances of these molecules at intermediate A_V. Thus, most of the columns of these molecules in molecular clouds are neither provided by the interior regions, where they are frozen out, nor by the surface, where they are photodissociated. Rather, most of the column arises from the intermediate zone.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 16 and 17 show the gas-phase elemental abundance of C and the gas-phase abundance ratio of elemental C/O, as a function of A_V, for two different cloud ages, t = 2 × 10⁶ years and t = 2 × 10⁷ years, respectively. In each figure, we plot the standard case, the standard case but with the density changed to n = 10³ and 10⁵ cm⁻³, and the standard case but with the incident FUV field changed to G₀ = 1 and 10³. For these results, we have used our time-dependent models, since the interesting effect of C/O exceeding unity occurs at high A_V, where the time-dependent models are more relevant and are more accurate in treating the carbon chemistry. At early times, t ∼ 2 × 10⁶ years, the C/O ratio generally stays below unity at high A_V. The exception is the high density (n = 10⁵ cm⁻³) case, but here the gas-phase elemental abundance of carbon is extremely low, which implies low abundances of any organic molecules produced in this carbon-rich environment. At late times, t ∼ 2 × 10⁷ years, we see that all cases except the high G₀ = 10³ case result in gas-phase C/O > 1 for A_V ≳ 8, due to the preferential freeze-out of water ice which robs the gas of elemental O. For those cases with C/O >1, the gas-phase elemental C abundance increases with increasing G₀ and decreasing n, due to increased thermal desorption rates of carbon molecular ices and decreased freeze-out rates, respectively.¹⁵ Figure 17 shows that for densities n ∼ 10³–10⁴ cm⁻³ and G₀ = 1–100, strong chemical effects might be observed deep in clouds, with high abundances of gas-phase hydrocarbons and other carbon molecules rather than CO. Unfortunately, these results can only be viewed as suggestive, because the gas-phase C depends on the desorption rate of CO, CH₄, and other carbon species frozen to grain surfaces, and these are uncertain. Nevertheless, it appears that regions with C/O >1 and high gas-phase elemental C abundance can appear deep in the cloud. This is consistent with the chemistry of GMC cores such as Orion A (e.g., Bergin et al. 1997) and regions in the Taurus Molecular Cloud (e.g., Pratap et al. 1997). In addition, in these deep layers significant amounts of CH₄ ice may form, and its release to the gas stimulates the carbon chemistry. In summary, the freeze-out of the oxygen into water ice can have profound effects on the carbon chemistry deep in the cloud and the ability of the cloud to produce carbon chain molecules or hydrocarbons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because our model suggests a peak in the H₂O and O₂ abundance at intermediate depths into a molecular cloud, one might expect clouds to be limb-brightened if viewed with beams smaller than the cloud. We make here two caveats to this prediction. First, GMCs are clumpy, and the clumpiness means that A_Vf will occur at different physical depths (i.e., radii) in the clouds, depending on the particular radial ray through the cloud and where it intersects clumps. This will have the effect of smearing the limb-brightening. In addition, although the 557 and 1113 GHz H₂O transitions are "effectively thin," they do suffer numerous scattering events since the optical depth in the lines is greater than unity. This scattering will reduce the limb-brightening effect for the H₂O transitions. The O₂ does not suffer scattering, so it may be more promising. However, the sensitivity of Herschel may require large O₂ columns >10¹⁶ cm⁻², and these may only be achieved for high G₀. In such cases, the peak O₂ abundance occurs deep in the cloud, A_V ∼ 8, and makes limb brightening only possible to detect in very high column but extended nearby clouds illuminated by high fields. Nevertheless, Herschel observations of O₂ in clouds may find regions of limb-brightening.

One still would like to map the clouds, however, even if limb brightening is hard to detect. The regions of high G₀ may be fairly limited in spatial extent, and our models predict higher H₂O and especially O₂ columns in these regions, which Herschel may spatially resolve. In addition, the H₂O and O₂ plateaus extend for ΔA_V ∼ 3. For gas densities of n ∼ 10⁴ cm⁻³, this corresponds to ∼0.2 pc, which at 500 pc corresponds to an angular size ∼15. Therefore, for a PDR slab viewed edge-on, Herschel may resolve the plateau region.

Even for clouds that are too small in angular extent to map, our model can be tested by correlation studies of the H₂O and O₂ intensities with the intensities of rotational transitions of other molecules known to trace either the surface or the volume of the cloud. For example, ¹³CO and C¹⁸O are thought to generally¹⁶ trace the volume of the cloud, whereas PDR species, such as CN, trace the "surface" of the cloud (Fuente et al. 1995; Rodriguez-Franco et al. 1998). We are investigating such observational correlations in Melnick et al. (2000).

The diffuse clouds that we model, based on those observed by SWAS toward W51 and W49 where H₂O 557 GHz is seen in absorption, are characterized by G₀ ∼ 1–10, n ∼ 100 cm⁻³, and A_V ∼ 1–5 (these are sometimes called "translucent clouds"; Snow & McCall 2006). Diffuse clouds are illuminated on all sides by the diffuse interstellar FUV field, whereas our models are those of 1D slabs illuminated from one side. To roughly construct a model of a diffuse cloud of total thickness A_Vt, we take the results from A_V = 0–0.5A_Vt and assume a mirror image for the rest of the cloud from 0.5A_Vt to A_Vt.

We compare a single model fit of a diffuse cloud with the 54 and 68 km s⁻¹ velocity components observed by SWAS toward the background submillimeter source W49 (Plume et al. 2004), as well as in the single absorption component toward W51 (Neufeld et al. 2002). For all three of these components, both OH and H₂O column densities have been measured, with N(H₂O) ∼10¹³–3 × 10¹³ cm⁻² and N(OH)/N(H₂O) ∼ 3–4. Our diffuse cloud model with n = 10² cm⁻³ and G₀ = 1 gives N(H₂O) = 3 × 10¹³ cm⁻² and N(OH)/N(H₂O) = 3.5 for a total cloud extinction of A_Vt = 4.4 magnitudes. We find that assuming a photodesorption yield $Y_{{\rm OH},w} \simeq 2 \times Y_{{\rm H_2O}}$ provides the best fit to the observed N(OH)/N(H₂O) ratios. The same model also provides predicted ¹²CO and atomic C columns and line intensities of the ¹²CO J = 1 − 0, J = 3 − 2, and [CI] 609 μm lines, all of which were observed in the W49 components. The models overpredict the columns and the intensities of the CO J = 1 − 0 and [CI] lines by a factor ∼2 that may be explained by dilution of the emission lines in the beam. The J = 3 − 2 line, however, is predicted to be a factor ∼2 smaller than is observed; a better fit for this line is found at higher densities, but makes the fit of all the other observed lines poorer. We note that Bensch et al. (2003) had similar problems in fitting a PDR model to the diffuse cloud CO and [CI] results.

Comparing our H₂O column to the model A_Vt, we find a column-averaged H₂O abundance of 3 × 10⁻⁹ with respect to H nuclei. The observational papers cited above quote H₂O/H₂ abundance ratios >10⁻⁷. These high abundances, however, were derived from the observed CO columns by multiplying the CO columns by 10⁴ to obtain H₂ columns. In our models, we find H₂ columns of ∼4 × 10²¹ cm⁻² and H₂O/H₂ abundance ratios of ∼10⁻⁸. The H₂/CO column-averaged abundance ratio in the model is 1.3 × 10⁵. In these translucent clouds, CO is photodissociated much more than H₂, which self-shields much more effectively. Considerable H₂ exists in regions where the gas phase carbon is mostly C⁺. These translucent clouds are examples of "dark H₂" which is not traced by CO.

4.2.1. A_V Thresholds for Ice Formation

4.2.2. B68

4.2.3. NGC 2024

4.2.4. O₂ Observed by Odin in ρ Oph

4.2.5. Upper Limits on O₂ Observed by Odin