Modeling Collisional Excitation of [O i] Fine Structure Line Emission from PDRs. I. Homogeneous Clouds

The energy from star formation put into a neighboring photon-dominated region (PDR) controls subsequent star formation in its surroundings, and it is thus critical to have reliable information on PDR structure and kinematics, as well as on the radiation from them, in order to understand how feedback from massive, young stars determines galactic evolution. Oxygen is the third most abundant element in the universe, and in a wide variety of astronomical environments it is found in neutral atomic form (O⁰). Fine structure line emission ([O I]) from atomic oxygen is the dominant coolant of warm PDRs and plays a key role in determining their structure and luminosity. The 63 μm line can be the most powerful far-infrared line in galaxies having relatively warm PDRs heated by massive star formation (Malhotra et al. 2001). Accurate analysis of [O I] emission has been hampered by the evident high optical depth of the 63 μm line and by the lack of a convenient and efficient code that can model radiative transfer in clouds with varying conditions along the line of sight.

The two lines at wavelengths 63.2 and 145.5 μm are widely used to analyze conditions in interstellar clouds and galaxies (e.g., Poglitsch et al. 1995; Oberst et al. 2011; Díaz-Santos et al. 2017). Almost all observations of these lines have been made with relatively low spectral resolution, so line profiles could not be determined. Several studies of regions of massive star formation carried with modest spectral resolution suggested [O I] absorption by foreground material (Poglitsch et al. 1996; Herrmann et al. 1997), but were not conclusive.

Observed line intensities and line profiles from some of the earliest detections (Stacey et al. 1983, 1993; Boreiko & Betz 1996) in giant molecular clouds suggested that the 63 μm line is often optically thick. This conclusion was consistent with straightforward single-component modeling that accompanied the observational results. Observational studies of regions of less-massive star formation and thus potentially with smaller gradients in conditions were surprising in terms of yielding extremely large values of line ratio I(145)μm/I(63)μm ≃ 0.4 (Caux et al. 1999; Liseau et al. 1999). These and the following papers cited reported fluxes (erg s⁻¹ cm⁻²) and ratios. We have converted these to intensities (erg s⁻¹ cm⁻² sr⁻¹) and calculated ratios assuming that beams are uniformly filled.

Their analyses indicated large optical depths in the 63 μm line, but the modeling with uniform physical conditions did not result in lower O⁰ column densities; Caux et al. (1999) concluded that N(O⁰)/N(CO) = 50, far greater than predicted by any chemical models. Saraceno et al. (1998) observed more than 70 pre-main-sequence objects with Infrared Space Observatory (ISO) and found higher I(145)μm/I(63)μm ratios than could readily be explained by models. These authors concluded that the explanation is colder gas associated with the source, rather than unrelated foreground material along the line of sight. Liseau et al. (2006) analyzed archival ISO data on many young stellar objects observed in both [O I] lines and suggested that tenuous foreground gas is absorbing the 63 μm line, making quantitative conclusions difficult. Nisini et al. (2015) observed jets in six young stellar objects, again finding relatively high I(145)μm/I(63)μm ratios. The jet emission lines are very broad (δv between 125 and 184 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$), so the suggestion of foreground absorption being responsible for high I(145)μm/I(63)μm ratios seems problematic. Their conclusion that H₂ molecule collisions are an order of magnitude less effective than H atoms in achieving a particular line ratio differs sharply from what we find with the more modern collision rate coefficients discussed in Section 2.

The upGREAT instrument on the Stratospheric Observatory for Infrared Astronomy (SOFIA) now offers the possibility of velocity-resolved spectroscopy of both [O I] lines, and recent observations (e.g., Leurini et al. 2015; Schneider et al. 2018; Mookerjea et al. 2019) have confirmed earlier suggestions (above references and Boreiko & Betz 1996) that optical depth effects are significant. [O I] spectra from SOFIA as well as those anticipated from the upcoming balloon mission GUSTO (https://www.nasa.gov/press-release/nasa-selects-mission-to-study-churning-chaos-in-our-milky-way-and-beyond) suggest that modeling tools that can predict line profiles as well as intensities will be very valuable. Anticipated new far-infrared spectroscopic observatories, including the ESA/JAXA M-Class Mission SPICA (Roelfsema et al. 2018) and the Origins Space Telescope proposed to the 2020 Decadal Survey (Leisawitz et al. 2018), will cover these spectral lines, albeit with a velocity resolution insufficient for detailed line profile studies of Galactic sources.

This paper is intended to take a step in this direction by presenting model results obtained with the Molpop–CEP statistical equilibrium/radiative transfer code (Elitzur & Asensio Ramos 2006; Asensio Ramos & Elitzur 2018). This code overcomes a major limitation of the commonly used large velocity gradient (LVG) model (e.g., Goldreich & Kwan 1974), which assumes that as a result of a large spatial velocity gradient, the radiative transfer is purely local, so that different portions of the spectral line are produced in different portions of the cloud. While satisfactory for many purposes and extremely computationally efficient, the LVG approximation is not able to model regions with very different physical conditions having the same velocity. Such regions are thus able to exchange photons, resulting in, for example, self-absorption features. Radiative transfer models that can handle interchange of photons more generally (e.g., Rybicki & Hummer 1991) are much more computationally demanding than LVG. The Molpop–Cep code achieves accuracy comparable to state-of-the-art accelerated Λ-iteration (ALI) codes while being orders of magnitude more rapid for cases with the many zones required to handle high optical depths (see discussion in Elitzur & Asensio Ramos 2006).

This paper is organized as follows. In Section 2 we discuss the fine-structure energy level structure of [O I] and recent results on collisions with molecular and atomic hydrogen, and in Section 2.2 we describe the behavior of the lines in the optically thin and optically thick limits. In Section 3 we present representative line profiles for a cloud with uniform conditions, and in Section 4 we give results for line intensities and line ratios as a function of H₂ density, kinetic temperature, and atomic oxygen column density.

The oxygen atom in its ground electronic state has three fine structure levels. As shown in Figure 1, they are inverted in energy, with the ${}^{3}{P}_{0}$ level having the highest energy and the ³P₂ the lowest. The two allowed transitions are also indicated; the ³P₀–³P₂ transition has a spontaneous decay rate approximately a factor of 10⁵ smaller and can be ignored for astronomical purposes. Information about the three fine structure transitions is given in Table 1. The collisional deexcitation rate coefficients are for a kinetic temperature of 100 K from Lique et al. (2018). These authors calculate the rate coefficients for H₂ in the J = 1 and J = 0 states individually, as well as the rates for collisions with H and He.

Zoom In Zoom Out Reset image size

The deexcitation rates for [O I] by collisions with atomic hydrogen are (as discussed by Lique et al. 2018) significantly smaller, by factors of 3 to 6, than those calculated previously by Abrahamsson et al. (2007). The Lique et al. (2018) deexcitation rate coefficients for collisions with H₂ are 2–4 times larger than those calculated by Monteiro & Flower (1987) for collisions with He producing the same transitions.

The deexcitation rates for collisions with H₂ are within a factor of 2 of those calculated by Jaquet et al. (1992) for deexcitation by ortho- and para-H₂. The values given here for collisions with H₂ are the average of the rates for the two spin states. The two rate coefficients for collisions to the ground state are comparable, and the rates for H₂ in J = 1 and J = 0 differ by less than 20%. The value for the ³P₀–³P₁ transition is approximately a factor of 50 smaller than the other two deexcitation rate coefficients, due to an approximate selection rule that forbids this transition (see discussion in Lique et al. 2018 with references to earlier treatments, including Monteiro & Flower 1987). This rate is a factor of ∼3 larger for H₂ in J = 1 than in J = 0, but the small magnitudes make using the average value of H₂ in J = 0 and J = 1 at all temperatures an acceptable approximation.

The major difference between collisions with H₂ and with H is the absence of the dramatically smaller ³P₀–³P₁ deexcitation rate coefficients for collisions with atoms, due to the presence of an indirect channel (Lique et al. 2018). Given the composition of the regions where the abundance of [O I] is expected to be appreciable, we focus on collisions with molecular hydrogen.

The rate equations for a three-level system presented and discussed in Goldsmith et al. (2015) can be applied to the [O I] fine structure lines, but with the caution that the ordering by total angular momentum is inverted for [O I] compared to that for [N ii]. The relevant expressions are given in the Appendix.

2.1. Relationships

Observations of the [O I] fine structure lines have been carried out by low-resolution spectrometers, such as PACS on Herschel (Poglitsch et al. 2010), and with high-spectral-resolution instruments that resolve line profiles (Risacher et al. 2016). While there is obviously significant additional information in the velocity-resolved spectra, current and future instrumentation, especially for extragalactic observations, will continue to employ the most sensitive spectrometers, which will not resolve the emission. In this paper, we give results for both types of observations, as well as some relationships for conversion between the two generally employed ways to present data.

Velocity-resolved spectroscopy has largely been carried out with single-mode diffraction-limited systems. The system is characterized by its effective area A_e, which applies to radiation coming from the direction of maximum sensitivity (the optical axis or boresight direction). The response to radiation coming from other directions is described by the normalized power pattern P_n, which is a function of the angle away from the boresight. Note that P_n has value unity for radiation incident from the boresight direction. We let such a system observe a source that has unpolarized specific intensity I_ν (units of, e.g., erg s⁻¹ cm⁻² sr⁻¹ Hz⁻¹), which is a function of angle ${\rm{\Omega }}^{\prime}$. Within a frequency interval δν at frequency ν, over which I is assumed to be constant, the power received in a single polarization with the telescope boresight at angle Ω is

$$\begin{eqnarray}&&{P}_{\mathrm{rec}}({\rm{\Omega }})=\displaystyle \frac{1}{2}{A}_{e}\delta \nu \int {P}_{n}({\rm{\Omega }}^{\prime} -{\rm{\Omega }}){I}_{\nu }({\rm{\Omega }}^{\prime} )d{\rm{\Omega }}^{\prime} .\end{eqnarray} \tag{ 1 }$$

If the source is uniform over the angular extent covered by P_n (or close to this limit), the expression simplifies to

$$\begin{eqnarray}&&{P}_{\mathrm{rec}}({\rm{\Omega }})=\displaystyle \frac{1}{2}{A}_{e}{I}_{\nu }\delta \nu \int {P}_{n}({\rm{\Omega }}^{\prime} )d{\rm{\Omega }}^{\prime} .\end{eqnarray} \tag{ 2 }$$

The integral in the above expression is the antenna solid angle, Ω_A. From the antenna theorem (e.g., Goldsmith 2002), for any single-mode diffraction-limited antenna, A_eΩ_A = λ². This gives us

$$\begin{eqnarray}&&{P}_{\mathrm{rec}}=\displaystyle \frac{{\lambda }^{2}}{2}{I}_{\nu }\delta \nu .\end{eqnarray} \tag{ 3 }$$

The power received from the antenna in a single-mode system defines the antenna temperature, T_A, through

$$\begin{eqnarray}&&{P}_{\mathrm{rec}}={{kT}}_{A}\delta \nu .\end{eqnarray} \tag{ 4 }$$

This gives the relationship for conversion between antenna temperature and specific intensity:

$$\begin{eqnarray}&&{T}_{A}=\displaystyle \frac{{\lambda }^{2}}{2k}{I}_{\nu }.\end{eqnarray} \tag{ 5 }$$

Here, T_A is a function of frequency, varying across the spectral line being observed, although this is generally not expressed explicitly in such expressions. Equation (5) allows conversion between antenna temperature and specific intensity for resolved spectral lines.

If we do not resolve the line, then what we observe is the antenna temperature integrated over frequency, which is related to the intensity $I=\int {I}_{\nu }d\nu$ through

$$\begin{eqnarray}&&\int {T}_{A}d\nu =\displaystyle \frac{{\lambda }^{2}}{2k}I.\end{eqnarray} \tag{ 6 }$$

Although the natural units for integrated antenna temperature are K Hz, it usually is expressed in units of K km s⁻¹, and we write

$$\begin{eqnarray}&&\int {T}_{A}(v){dv}\,({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1})=\displaystyle \frac{{\lambda }^{3}}{2\times {10}^{5}k}I.\end{eqnarray} \tag{ 7 }$$

Table 2 gives values for the constants relating antenna temperature and intensity for the two [O I] fine structure transitions. An approximate value for the antenna temperature can be obtained from the intensity by assuming a simple line shape and line width Δv. Then, dividing the value of $\int {T}_{A}{dv}$ obtained using the fourth column of Table 2 by the assumed Δv yields the peak antenna temperature.

2.2. Limiting Cases

It is valuable to consider some limiting cases of excitation of the [O I] fine structure lines. For this purpose, we use analytic expressions that give insight into the behavior of the system. Equivalent expressions have been previously developed for three level systems in articles (e.g., Liseau et al. 2006) and books (Draine 2011), but we start here from basic principles for completeness and clarity.

Here, and in what follows, we use the superscript 63 to denote the ³P₁–³P₂ transition having a wavelength ≃63 μm, and 145 to denote the higher-lying ³P₀–³P₁ transition with λ ≃ 145 μm.

We denote the collision rate coefficient from level i to level j as R_ij, which has units cm⁻³ s⁻¹. The rate coefficients are obtained by integrating the cross sections over the velocity distribution of the colliding particles. The R_ij values here are denoted γ_ij by Spitzer (1978) and q_ij by Osterbrock (1989). For collisions with species k (here hydrogen atoms or hydrogen molecules), the collision rate is C_ij = R_ijn(k), having units s⁻¹. The rate of collisional transitions per unit volume is equal to n_i(O)n(k)R_ij = n_i(O)C_ij.

2.2.1. Optically Thin Subthermal Emission

Subthermal excitation occurs when the downward collision rates are much smaller than competing radiative rates. In this limit, the collisional processes that enter into the rate equations for the level populations (density of oxygen in level i, n_i(O)) are exclusively upward (excitation) collisions, and we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{3}}{{n}_{2}}=\displaystyle \frac{{A}_{21}}{{A}_{32}}\displaystyle \frac{{C}_{13}}{{C}_{12}+{C}_{13}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{2}}{{n}_{1}}=\displaystyle \frac{{C}_{12}+{C}_{13}}{{A}_{21}},\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&{n}_{1}={n}_{t}{\left[1+\displaystyle \frac{{C}_{12}+{C}_{13}}{{A}_{21}}+\displaystyle \frac{{C}_{13}}{{A}_{32}}\right]}^{-1},\end{eqnarray} \tag{ 10 }$$

where n_t is the total volume density of O⁰ atoms. These results depend on the relative magnitudes of the different collision rate coefficients. Note that C₂₃ does not appear, and that the ratio n₃/n₂ is independent of density in the subthermal limit, while the ratio n₂/n₁ is proportional to the collision rate and thus the density. From Table 1, the ratio A₂₁/A₃₂ is equal to 5.1, and with the statistical weights of level i, g_i = 2i + 1 for i = 1, 2, and 3, the upward collision rates (for kinetic temperature 100 K) are C₁₂ = 1.1 × 10⁻¹¹ n(H₂) s⁻¹ and C₁₃ = 1.0 × 10⁻¹² n(H₂) s⁻¹. We then find the population ratios n₂/n₁ = 1.4 × 10⁻⁷ n(H₂) and ${\text{}}{n}_{3}/{n}_{2}=0.43$.

Table 1.  [O I] Fine Structure Transitions and Collisional Parameters

Transition	Frequencya	Wavelength	${E}_{{\rm{u}}}$/k	Aulb	Rul(H)c	Rul(H2)c
	(GHz)	(μm)	(K)	(s−1)	(10−10 cm3 s−1)	(10−10 cm3 s−1)
3P0–3P1	2060.069	145.53	326.6	1.7 × 10−5	0.84	.0291
3P1–3P2	4744.777	63.18	227.7	8.7 × 10−5	1.12	1.74
3P0–3P2	6804.847	44.06	326.6	1.4 × 10−10	0.76	1.36

Notes.

^aFrom Zink et al. (1991); these values supersede those of Saykally & Evenson (1979). ^bFrom Fischer & Saha (1983). There are slight differences among different calculations and references; see Baluja & Zeippen (1988). ^cAt kinetic temperature 100 K.

Download table as:  ASCIITypeset image

The excitation temperature describes the relative populations per unit statistical weight and for transition from level u to level l is defined as

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}^{{ul}}=\displaystyle \frac{{T}_{{ul}}^{* }}{\mathrm{ln}({n}_{l}{g}_{u}/{n}_{u}{g}_{l})},\end{eqnarray} \tag{ 11 }$$

where ${T}_{{ul}}^{* }$ = ΔE_ul/k is the equivalent temperature of the transition with ΔE the energy difference between the upper and lower levels. This yields the subthermal limits for the excitation temperatures (for collisions with H₂ at a kinetic temperature of 100 K) T_ex¹⁴⁵ = −387 K, and for a density of 10³ cm⁻³, T_ex⁶³ = 27 K.

The results from Molpop–CEP with an O⁰ column density sufficiently small that optical depths are far less than unity are shown in Figure 2. The 63 μm transition is well behaved, with T_ex⁶³ rising smoothly toward the 100 K kinetic temperature. Here, T_ex⁶³ for n(H₂) = 10³ cm⁻³ is again 27 K. The figure includes results for collisions with atomic hydrogen (Lique et al. 2018). The O–H rates being a few times smaller than those for O–H₂ collisions results in slightly lower values of T_ex at intermediate densities.

Zoom In Zoom Out Reset image size

From the full calculation, at densities n(H) ≤ 10² cm⁻³, T_ex¹⁴⁵ = −333 K, close to the approximate value of −387 K given above and representing a very small difference in the level populations. The population inversion of the 145 μm transition (noted previously with earlier collision rate coefficients by, e.g., Liseau et al. 2006 and with more recent coefficients by Lique et al. 2018) is evident for densities n(H₂) ≤ 10⁵ cm⁻³. For this range of densities, the rapid depopulation of the lower (³P₁) level by spontaneous decay to the ground state overwhelms the combination of radiative decay from the upper (³P₀) level and collisional population of the lower level by collisions from the ground (³P₂) state.

The presence of the population inversion is due to the relative magnitude of the rates of populating the levels by collisions compared to depopulation by spontaneous emission. The relative collision rate coefficients depend only on the kinetic temperature. As discussed above, at 100 K the ratio n₃/n₂ = 0.43, which is larger than the ratio of the statistical weights g₃/g₂ = 0.33. In order to eliminate the population inversion, we require that C₁₂/C₁₃ (equal to 10.4 at 100 K) be greater than or equal to 14.6. At a kinetic temperature of 50 K, for example, this ratio is ≃28, because the rate from the ground state to the second excited state is reduced more by the low kinetic temperature than is the rate from the ground state to the first excited state. The result is indeed no population inversion, but rather a modestly superthermal excitation with T_ex¹⁴⁵ ≃ 150 K.

The inversion potentially provides amplification, but the column density of atomic oxygen in plausibly sized regions of this low density and the suppression of the inversion by the buildup of intensity at the frequency of the transition together prevent the development of any anomalously strong line intensities. For densities greater than 10⁵ cm⁻³, collisional deexcitation becomes significant and Equation (8) is no longer valid.

The intensity of optically thin emission in a transition from level u to level l is given by

$$\begin{eqnarray}&&{I}^{{ul}}=\displaystyle \frac{{A}_{{ul}}{N}_{u}h\nu }{4\pi }.\end{eqnarray} \tag{ 12 }$$

The ratio of intensities of two transitions can yield valuable information. Assuming that the ratio of the column densities is equal to that of the volume densities, we find for optically thin emission from [O I] that

$$\begin{eqnarray}&&\displaystyle \frac{{I}^{145}}{{I}^{63}}=\displaystyle \frac{{A}_{32}}{{A}_{21}}\displaystyle \frac{{\nu }_{32}}{{\nu }_{21}}\displaystyle \frac{{n}_{3}}{{n}_{2}}=0.085\displaystyle \frac{{n}_{3}}{{n}_{2}}.\end{eqnarray} \tag{ 13 }$$

In the subthermal limit, we use Equation (8) and obtain

$$\begin{eqnarray}&&\displaystyle \frac{{I}^{145}}{{I}^{63}}=\displaystyle \frac{{\nu }_{32}}{{\nu }_{21}}\displaystyle \frac{{C}_{13}}{{C}_{12}+{C}_{13}}=\displaystyle \frac{0.434}{1+\tfrac{{C}_{12}}{{C}_{13}}}.\end{eqnarray} \tag{ 14 }$$

The line intensity ratio I¹⁴⁵/I⁶³ for optically thin emission increases as the kinetic temperature increases in the subthermal (low excitation rate) limit, growing by a factor of ≃5 as the kinetic temperature goes from 50 to 400 K. This temperature dependence is a reflection of the greater energy requirement for the primary excitation path for the 145 μm line, C₁₃, compared to that for the 63 μm line, C₁₂. If the volume density can be constrained to ≤10⁵ cm⁻³ and the atomic oxygen column density to ≤10¹⁷ cm⁻², the line ratio can be a useful probe of the temperature in the PDR.

For a kinetic temperature of 100 K, the ratio of the intensities is

$$\begin{eqnarray}&&\displaystyle \frac{{I}^{145}}{{I}^{63}}{| }_{\mathrm{low} \mbox{-} n100{\rm{K}}}=0.038.\end{eqnarray} \tag{ 15 }$$

The intensity can be converted to peak or integrated intensity through Equations (5) to (7). The antenna temperature ratio of the two [O I] transitions, if both are optically thin and have equal line widths, can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{A}^{145}}{{T}_{A}^{63}}={\left(\displaystyle \frac{{\lambda }_{32}}{{\lambda }_{21}}\right)}^{2}\displaystyle \frac{{A}_{32}}{{A}_{21}}\displaystyle \frac{{n}_{3}}{{n}_{2}}.\end{eqnarray} \tag{ 16 }$$

In the subthermal limit, we find from Equations (8) and (16) that

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{A}^{145}}{{T}_{A}^{63}}=5.31\displaystyle \frac{{C}_{13}}{{C}_{13}+{C}_{12}}=\displaystyle \frac{5.31}{1+\tfrac{{C}_{12}}{{C}_{13}}},\end{eqnarray} \tag{ 17 }$$

and for low H₂ densities at a kinetic temperature of 100 K, we find

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{A}^{145}}{{T}_{A}^{63}}{| }_{\mathrm{low} \mbox{-} n100{\rm{K}}}=0.47.\end{eqnarray} \tag{ 18 }$$

The primary dependence on kinetic temperature is through the difference in upper state energies as included in Equation (17), with weaker dependence from the temperature dependence of the collision rates themselves.

2.3. Optical Depth

The optical depth of a spectral line is proportional to its line profile function ϕ(ν) describing its frequency dependence, normalized such that $\int \phi (\nu )d\nu =1$. The optical depth τ is given by

$$\begin{eqnarray}&&\tau (\nu )=\displaystyle \frac{h\nu }{c}\phi (\nu )({N}_{l}{B}_{{lu}}-{N}_{u}{B}_{{ul}}),\end{eqnarray} \tag{ 19 }$$

where N_l and N_u are the lower- and upper-level column densities, respectively, and B_lu and B_ul are the upward and downward stimulated emission coefficients. Assuming a Gaussian line profile with dispersion σ and FWHM line width Δν = 1.67σ, we find that the line profile function at the line center is $\phi ({\nu }_{0})=1/(\sqrt{\pi }\sigma )=0.94/{\rm{\Delta }}\nu$. Using the line width expressed in $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, we obtain

$$\begin{eqnarray}&&\tau ({\nu }_{0})=\displaystyle \frac{0.94{{hB}}_{{lu}}{N}_{l}}{{10}^{5}{\rm{\Delta }}v(\mathrm{km}\ {{\rm{s}}}^{-1})}(1-{e}^{-{T}^{* }/{T}_{\mathrm{ex}}}).\end{eqnarray} \tag{ 20 }$$

Using the relationship between upward and downward stimulated emission coefficients and the relationship between the B and A coefficients, we find

$$\begin{eqnarray}&&\tau ({\nu }_{0})=\displaystyle \frac{0.94{A}_{{ul}}{\lambda }^{3}}{8\pi {10}^{5}{\rm{\Delta }}v(\mathrm{km}\ {{\rm{s}}}^{-1})}\displaystyle \frac{{g}_{u}}{{g}_{l}}{f}_{l}N(1-{e}^{-{T}^{* }/{T}_{\mathrm{ex}}}),\end{eqnarray} \tag{ 21 }$$

where N is the total column density of the species and f_l is the fraction in the lower level of the transition.

We define a maximum line center optical depth, τ₀, by taking the FHWM line width equal to 1 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, a lower-level column density of 1 cm⁻², and highly subthermal excitation so that ${T}^{* }/{T}_{\mathrm{ex}}$ ≫ 1, allowing us to neglect the exponential term. This yields

$$\begin{eqnarray}&&{\tau }_{0}=\displaystyle \frac{0.94{A}_{{ul}}{\lambda }^{3}}{8\pi {10}^{5}}\displaystyle \frac{{g}_{u}}{{g}_{l}}\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&\tau ({\nu }_{0})={\tau }_{0}\displaystyle \frac{{f}_{l}N({\mathrm{cm}}^{-2})}{{\rm{\Delta }}v(\mathrm{km}\ {{\rm{s}}}^{-1})}.\end{eqnarray} \tag{ 23 }$$

Table 3 gives values of τ₀ and f_l for local thermodynamic equilibrium (LTE) at different temperatures. Even with the assumption of LTE, the fractional population of the lower level of the 63 μm line (³P₂ level) is not far from unity, and that for the lower level of the 145 μm line (³P₁ level) is between 0.02 and 0.07.

PDRs are material subjected to ultraviolet and visible radiation that cannot ionize hydrogen, but deposits energy that both heats the gas and dust and controls the chemical composition of these regions (Tielens & Hollenbach 1985). Figure 3 shows aspects of a PDR relevant to [O I] emission, generated using the Meudon PDR code (Le Petit et al. 2006). In a PDR exposed to a radiation field enhanced by a factor of 10⁴ relative to the standard interstellar radiation field (which is easily attained in regions of massive star formation) and density of 10⁴ cm⁻³, oxygen is essentially all in atomic form for visual extinctions A_v ≤ 3.5 and is still over 30% atomic for A_v = 10 mag.

Zoom In Zoom Out Reset image size

The transition from atomic to molecular hydrogen with the parameters of the present model occurs at A_v ≃ 0.9 mag. At a density of 10⁴ cm⁻³, the kinetic temperature is somewhat over 500 K for A_v ≤ 0.4 mag, drops to 200 K at A_v = 1.6 mag, and continues to drop beyond that point. Even in the higher temperature portion of this PDR, the fraction of oxygen in the ³P₁ level is only a few percent, and far less than this when the temperature drops below 100 K. The total column density of atomic oxygen in the PDR is on the order of 3 × 10¹⁸ cm⁻². We thus can take f_l = 1 for the 63 μm line and f_l = 0.01 for the 145 μm line. From Table 3 and Equation (23) we find that, for a line width of 5 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and the above parameters, τ¹⁴⁵ (ν₀) = 0.04 and τ⁶³ (ν₀) = 3. It is readily apparent that the 63 μm line can be significantly optically thick, while the 145 μm line will almost certainly be optically thin.

2.4. [O i] Absorption

2.5. Transition to the Optically Thick Regime

The radiative trapping resulting from a transition being optically thick can be thought of as reducing the spontaneous emission rate, thus bringing the excitation closer to LTE. This is often relevant for the 63 μm transition, which is subthermally excited for a wide range of densities up to n(H₂) ≃ 10⁶ cm⁻³ (see discussion in Section 2.6). Figure 4 shows this behavior for atomic oxygen in a region with n(H₂) = 10³ K. As shown by the red curve in the lower panel, the optical depth of the 63 μm transition increases almost linearly with O^o column density and reaches a value of 10 for N(O^o) = 4 × 10¹⁸ cm⁻². At this point, the excitation temperature indicated by the red curve in the upper panel has increased appreciably above the value ≃27 K produced by collisional excitation alone.

Zoom In Zoom Out Reset image size

The upper panel of Figure 4 shows an additional radiative effect. In this implementation of the Molpop–CEP program, there is no systematic velocity shift through the slab of material being modeled. Thus, the oxygen atoms in the center of the slab are exposed to a stronger radiation field at a frequency of 4475 GHz than are the atoms at its edge. The excitation temperature at the center of the slab (denoted C) is significantly greater than that at the edge (denoted E). This difference increases throughout the range of optical depths included here (up to ≃200), but for sufficient oxygen column density (greater than included here), the transition becomes thermalized throughout.

The situation for the 145 μm transition is more complex, as the adopted density of 10³ cm⁻³ produces a population inversion (negative excitation temperature) in the optically thin limit. As the column density increases, so does τ¹⁴⁵, but its magnitude is very small, ≤10⁻⁵. When N(O^o) reaches 10¹⁷ cm⁻², the excitation temperature switches to being positive but very large, indicating essentially equal populations in the upper and lower levels. For higher column densities, T_ex¹⁴⁵ drops, indicating a reduction in the ratio n₃/n₂. This is not due to the trapping in the 145 μm line, for which the optical depth is very small, ≃10⁻⁵–10⁻². Rather, the trapping in the 63 μm line effectively reduces its spontaneous decay rate A₂₁, which in turn reduces the the ratio n₃/n₂ given by Equation (8), and hence T_ex¹⁴⁵.

We can gain a feeling for this effect by assuming that the reduction in the A coefficient is given by the escape probability β = (1 − exp − τ)/τ ≃ 1/τ for τ ≫ 1. Then for τ⁶³ = 10, βA₂₁/A₃₂ is reduced to 0.51 from its optically thin value of 5.1. The result from Equations (8) and (11) is that T_ex¹⁴⁵ ≃ 50 K. This is intermediate between the C and E values in Figure 4. Note that the value of T_ex¹⁴⁵ drops to values below 30 K as a consequence of this effect. Only for oxygen column densities greater than ∼10²⁰ cm⁻³ does τ¹⁴⁵ become sufficiently large that the trapping in this transition begins to drive the excitation temperature up toward the kinetic temperature. The conclusion is that for extended low-density regions with atomic oxygen column densities between 10¹⁸ cm⁻² and as high as a few ×10²⁰ cm⁻², the 145 μm transition can have an excitation temperature below that of the 63 μm transition, with both lines being optically thick.

2.6. Critical Density

The critical density is defined as the density of colliding partners that makes the downward collision rate equal to the spontaneous decay rate, so

$$\begin{eqnarray}&&{R}_{{ul}}{n}_{c}={A}_{{ul}}.\end{eqnarray} \tag{ 24 }$$

As discussed by Goldsmith et al. (2015), in a two-level system, this makes the ratio of upper- to lower-level populations equal to one-half of the value in LTE at the kinetic temperature. Using the definition of excitation temperature (Equation (11)), we define the critical excitation temperature, T_ex^c, as that resulting from a density equal to the critical density. In general,

$$\begin{eqnarray}&&\displaystyle \frac{\tfrac{{n}_{u}}{{n}_{l}}}{{\tfrac{{n}_{u}}{{n}_{l}}}_{\mathrm{LTE}}}={e}^{-{T}^{* }/{T}_{k}+{T}^{* }/{T}_{\mathrm{ex}}},\end{eqnarray} \tag{ 25 }$$

and from the definition of the critical density, we obtain

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}^{c}=\displaystyle \frac{{T}^{* }}{\tfrac{{T}^{* }}{{T}_{k}}+\mathrm{ln}2}.\end{eqnarray} \tag{ 26 }$$

From Table 1, the two-level critical density for the 63 μm [O I] line is n_c = 5.0 × 10⁵ cm⁻³, and from the above, T_ex^c = 77 K. From Figure 2, we see that the full calculation gives very close to this same value. The two-level and full calculations agree, due to the relatively low fractional population of the upper ³P₀ level and the unimportance of collisions involving that level in determining the properties of the lower ³P₁–³P₂ 63 μm transition.

The situation for the 145 μm line is more complex because there is a population inversion at low densities, so the excitation temperature approaches the kinetic temperature from above as the hydrogen density increases and we reach thermalization. For collisions with H₂, the unusually low collisional deexcitation rate coefficient (Table 1) gives a very large value of n_c = 5.8 × 10⁶ cm⁻³. At this density, the excitation temperature of the 145 μm line is only ≃10 K greater than the kinetic temperature. Thus, taking the critical density as defined by Equation (24) gives an unrealistically high value for the H₂ density required to bring the excitation temperature close to the kinetic temperature, as collisional pathways other than ³P₀–³P₁ dominate.

For collisions with atomic hydrogen, the ³P₀–³P₁ deexcitation rate coefficient has a more typical value, and n_c(H⁰) = 2 × 10⁵ cm⁻³. From the full calculation, T_ex^c is several times the kinetic temperature at this density, but is dropping rapidly toward the kinetic temperature as the density increases and we approach LTE. For both collision partners, the critical density for the 145 μm line has limited usefulness. Table 4 gives the values of the critical densities for the two [O I] fine structure transitions and collisions with H and H₂. The values for collisions with H^o are a factor of ≃2 larger than those given by Fernández-Ontiveros et al. (2016).

2.7. Optically Thin Thermalized Emission

We can use Equation (13) with the LTE value of the density ratio

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{3}}{{n}_{2}}=\displaystyle \frac{{g}_{3}}{{g}_{2}}{e}^{-{T}^{* 32}/{T}_{k}}\end{eqnarray} \tag{ 27 }$$

to obtain

$$\begin{eqnarray}&&\displaystyle \frac{{I}^{145}}{{I}^{63}}{| }_{\mathrm{LTE}}=0.0283{e}^{-98.8/{T}_{k}},\end{eqnarray} \tag{ 28 }$$

which gives

$$\begin{eqnarray}&&\displaystyle \frac{{I}^{145}}{{I}^{63}}{| }_{\mathrm{LTE}100{\rm{K}}}=0.011.\end{eqnarray} \tag{ 29 }$$

Figure 5 shows the variation of line intensity ratio as a function of molecular hydrogen density for four values of the kinetic temperature. The behavior for these is largely similar, with I¹⁴⁵/I⁶³ dropping by approximately a factor of 4 as the density increases and the excitation changes from subthermal to thermalized. The temperature dependence in the subthermal limit is given in Equation (14), somewhat different than that given in the thermalized limit by Equation (28). In the latter case, the ratio increases as the kinetic temperature increases, due to the higher energy of the upper level of the 145 μm transition, but is independent of the temperature dependence of the collision rates, which do enter for subthermal excitation.

Zoom In Zoom Out Reset image size

For 400 K kinetic temperature, we show in Figure 5 the results for collisions with atomic hydrogen as well as molecular hydrogen. Only a small difference can be discerned for low densities, at which the excitation is subthermal, while for higher densities, the two collision partners yield an identical ratio as the transitions are thermalized. This is quite different from the conclusion of Nisini et al. (2015), who found an order of magnitude or more difference in density to obtain a given intensity ratio when considering collisions with H₂ and H^o. This is presumably a result of different collision rate coefficients; those employed by Nisini et al. (2015) were not specified.

2.7.1. Optically Thick Thermalized Emission

In this limit, the antenna temperature is just that of a blackbody at the kinetic temperature T_k:

$$\begin{eqnarray}&&{T}_{A}=\displaystyle \frac{{T}^{* }}{{e}^{{T}^{* }/{T}_{k}}-1}.\end{eqnarray} \tag{ 30 }$$

The ratio of the antenna temperatures of the two [O I] transitions in this limit is

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{A}^{145}}{{T}_{A}^{63}}=\displaystyle \frac{{T}^{* 145}}{{T}^{* 63}}\displaystyle \frac{{e}^{{T}^{* 63}/{T}_{k}}-1}{{e}^{{T}^{* 145}/{T}_{k}}-1}=0.434\displaystyle \frac{{e}^{227.7/{T}_{k}}-1}{{e}^{98.8/{T}_{k}}-1}.\end{eqnarray} \tag{ 31 }$$

For a kinetic temperature of 100 K, T_A¹⁴⁵ = 58.6 K while T_A⁶³ = 26.0 K, and their ratio is

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{A}^{145}}{{T}_{A}^{63}}{| }_{\mathrm{thick}100{\rm{K}}}=2.25.\end{eqnarray} \tag{ 32 }$$

This high value is a reflection of the greater departure from the Rayleigh–Jeans limit of the shorter wavelength transition. In consequence, the ratio approaches unity for T_k ≥ 250 K, but can become very large for kinetic temperatures appreciably less than 150 K.

The profile of a spectral line can provide information about the line-of-sight variations in physical conditions if there are large-scale velocity gradients through the region of interest. The optical depth has a strong influence on what we can discern in terms of source structure. Optically thin emission plus velocity shifts encode the information about physical conditions into the velocity dependence of the signal. Optically thick emission without large-scale velocity shifts hides from us more distant portions of gas along the line of sight, but even if the emission is optically thick, we can gain some information, such as about the temperature. Optically thick emission can also function as a background source against which lower-excitation foreground gas can be seen in absorption.

Modeling the emergent line spectrum requires following the emission and absorption of photons of different frequencies as they propagate through the cloud. The Molpop–CEP program carries this out, without making any assumptions about the magnitude of velocity gradients (as does the large velocity gradient or LVG model). Photons of a given frequency can be absorbed or can stimulate emission, or can be added to by spontaneous emission, from any point in the cloud. To calculate accurately the line profile, we do have to input the physical conditions as a function of position in the cloud, which is here taken as a slab, with variations only along the line-of-sight coordinate. Doing this can be complex, and there is a real issue in terms of uniqueness in terms of inversion of the observed line profile to determine the source structure. The enormous information content of spectrally resolved observations nevertheless makes this a critical tool for really determining the physical parameters and their variation within a source.

The present paper does not attempt to build a detailed source model (as has been done by, e.g., Schneider et al. 2018). Rather, we give a few examples that indicate the variety of behaviors that can be expected even for a uniform slab model of an [O I] source, and to connect them to the discussion of general behavior given above. We present the results for line profiles in units of antenna temperature, but these can be easily converted to specific intensity through Equation (5) and the values given in Table 2.

Table 2.  Conversion Factors

Transition	ΔE/k	TA/Iν	$\int {T}_{A}{dv}/I$
	(K)	(K/erg s−1 cm−2 sr−1 Hz−1)	(K km s−1/erg s−1 cm−2 sr−1)
145	98.8	7.67 × 1011	1.12 × 105
63	227.7	1.45 × 1011	9.13 × 103

Download table as:  ASCIITypeset image

When both transitions are optically thin, the excitation temperature is purely a result of collisional excitation and is uniform throughout the slab. However, when the atomic oxygen column density increases so that the 63 μm line becomes optically thick, radiative trapping starts to affect the level populations and the excitation temperature. Since the cloud modeled here has no large-scale velocity gradient, this first occurs at the center of the cloud, where the effect of the optical depth is greatest.

Figure 6 shows the behavior of the excitation temperature as a function of position within the slab having N(O^o) = 10¹⁹ cm⁻² and FWHM line width  1.67 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, resulting in optical depths of 25–30 for the 63 μm transition and 0.2–1.5 for the 145 μm transition. We thus see significant variation in the excitation temperatures as a function of position in the slab for the two lower volume densities. The 63 μm line with no trapping has T⁶³_ex = 27 K (see Figure 2), but as seen in Figure 6, this has been increased to 50 K in the center of the slab, but drops to 39 K at the slab edge. Similar behavior occurs for n(H₂) = 10⁵ cm⁻³, but starting from a purely collisionally produced excitation temperature of 56 K (lower than that seen even in the outermost zone included in Figure 6 due to the radiation from the interior of the cloud).

Zoom In Zoom Out Reset image size

The 145 μm line is inverted (T¹⁴⁵_ex = −339 K) by collisions for n(H₂) = 10³ cm⁻³ and superthermal (T¹⁴⁵_ex = 1550 K) for n(H₂) = 10⁵ cm⁻³ (see Figure 2). The optical depth, τ⁶³ = 30, results in the relatively low excitation temperature T¹⁴⁵_ex = 30 K in the slab center (see discussion in Section 2.5), which rises as one moves toward the edges of the slab. For lower N(O^o), T¹⁴⁵_ex rises to above the kinetic temperature and becomes negative at the boundary of the slab. For n(H₂) = 10⁷ cm⁻³, both transitions are nearly thermalized by collisions alone, which is reflected in the minimal variation in excitation temperature throughout the slab.

Figure 7 shows the line profiles for the two [O I] fine structure lines for four values of the atomic oxygen column density. The kinetic temperature is 100 K, the H₂ density is 10⁷ cm⁻³, and the FWHM line width is 1.67 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in all cases. For the lowest column density, both lines are optically thin (τ⁶³ = 0.22), and we see nearly Gaussian line profiles. The ratio of peak antenna temperature is ${T}_{A}^{145}/{T}_{A}^{63}=0.25$, very close to the ratio from Equations (16) and (27) for quasi-thermalized excitation.

Zoom In Zoom Out Reset image size

As the column density increases, we see that the 63 μm line becomes flat-topped due to saturation for N(O^o) = 10¹⁹ cm⁻². For this and higher column densities, the peak antenna temperature of the 145 μm line is greater than that of the 63 μm line, due to the departure from the Rayleigh–Jeans limit discussed in Section 2.7.1. While the largest column density N(O^o) = 10²⁰ cm⁻² is not widespread, it can be found in massive PDRs (Goicoechea et al. 2015), and it illustrates the limit of a uniform cloud with both transitions thermalized and optically thick.

Figure 8 illustrates the more complex line profiles that can result when the density (and hence the collision rate) alone is insufficient to thermalize the transitions. In this case, with n(H₂) = 10⁴ cm⁻³ and T_k = 100 K, both transitions are, in the absence of trapping, subthermally excited. When the column density reaches 10¹⁸ cm⁻², the fractional population of the ³P₂ state, the lower level of the 63 μm transition, is larger than it is in LTE at T_k = 250 K, so the optical depth is larger, producing the slightly flat-topped appearance.

Zoom In Zoom Out Reset image size

For N(O^o) = 10¹⁹ cm⁻², we see for the first time clear evidence of the importance of radiative trapping and the effect of coupling the radiation from different regions within the slab. Under these conditions, the excitation temperature of the 63 μm line is 79 K in the center of the two-sided slab, but drops to ≤60 K in the outer 5% of the slab. Since the total optical depth through the slab of the line is 28, the outer 5% corresponds to optical depth ≃1. Thus, the antenna temperature at the line center reflects the lower value resulting from purely collisional excitation, while moving into the line wings, one sees deeper into the cloud where the radiative trapping increases the excitation temperature. There is sufficient opacity in the 63 μm line to result in higher antenna temperatures at velocities offset by ±1.5 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ than at line center.

The 145 μm line has a total optical depth of 0.9, but the radiative trapping here is having the effect, discussed earlier, of suppressing the population inversion produced by collisions. The result is that the excitation temperature T¹⁴⁵_ex is 44 K at the center of the slab, but is 50–60 K at the edge of the slab where the trapping is less. For N(O^o) = 10¹⁹ cm⁻², the line profile is close to Gaussian as the optical depth is modest.

For the large column density N(O⁰) = 10²⁰ cm⁻², both transitions are optically thick. The 63 μm line profile is similar to that for a column density an order of magnitude lower, but the central depression is wider because of the larger (τ⁶³ = 260) optical depth. The 145 μm line profile in this situation is similar to that of the shorter-wavelength line because the excitation temperature is 75 K in the center of the slab but only 56 K at the slab edge, with a total optical depth of 16 resulting again in a centrally depressed profile. The behavior of T¹⁴⁵_ex as a function of column density is shown in Figure 4; for n(H₂) = 10⁴ cm⁻³, its behavior is similar to that for n(H₂) = 10³ cm⁻³ illustrated, with the excitation temperature being driven to very low values by trapping of the 63 μm line, and with T¹⁴⁵_ex increasing toward the center of the slab, where trapping enhances the radiation field at the frequency of the 145 μm transition.

The results from models of uniform slabs depend on numerous parameters and can be presented in many ways. The overwhelming majority of data on [O I] fine structure lines to date has been obtained with systems having relatively low spectral resolution that do not resolve the line profiles. Such systems continue to be available (FORCAST on SOFIA) and will possibly be a major source of data in the future (SPICA). With this in mind, we first present some results for line intensities, even though the output from Molpop–CEP always provides the spectral line profiles.

4.1. Line Intensities

We first show results in terms of line intensities for four kinetic temperatures, with the intensity of the 63 μm line and the ratio of the intensity of the 145 μm line to that of the 63 μm line as a function of O⁰ column density, for specified values of the volume density, n(H₂).

The intensity of the 63 μm line (lower panels of Figures 9–12) increases monotonically with O⁰ column density and volume density. This line is quite optically thick (τ ≥ 3) for N(O⁰) ≥ 10¹⁸ cm⁻² (Table 3) and densities below the critical density (Table 4). The emitted intensity continues to increase even though the line is optically thick because collisional deexcitation is unimportant compared to spontaneous decay, and the emission is effectively optically thin (e.g., Goldsmith et al. 2012). At volume densities for which collisional deexcitation is important and the excitation temperature is comparable to the equivalent temperature of the transition, the line intensity saturates for N(O⁰) ≥ 10¹⁸ cm⁻², at which point the emission is thermalized and optically thick.

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

Table 3.  Maximum Line Center Optical Depth and Lower-level Fractional Populations for [O I] Fine Structure Transitions

Transition	τ0	fl
		T(K)
		50	100	250	400
145	6.52 × 10−18	6.4 × 10−3	5.8 × 10−2	1.9 × 10−2	2.4 × 10−1
63	4.92 × 10−18	9.9 × 10−1	9.4 × 10−1	7.7 × 10−1	7.0 × 10−1

Download table as:  ASCIITypeset image

Table 4.  Critical Densities for [O I] Fine Structure Transitions

Transition	nc(H2)	nc(H)
	(cm−3)	(cm−3)
145	5.8 × 106	2.0 × 105
63	5.0 × 105	7.8 × 105

Download table as:  ASCIITypeset image

The behavior of the line intensity ratio, shown in the upper panels of Figures 9 through 12, was discussed in Section 2.2. The temperature dependence of I¹⁴⁵/I⁶³ is similar in the thermalized and subthermal limits. We here see that the dependence on atomic oxygen column density is more complex. The density dependence in the high-column-density limit is inverted compared to that in the low-column-density limit. This is because at low volume densities, the ratio is "anomalously" high because of the peculiar combination of spontaneous emission rates and collision rate coefficients discussed in Section 2.2.1. As the collision rate increases, this ratio is lowered to that characteristic of LTE.

At high column densities, radiation trapping in the optically thick 63 μm line reduces the population of the uppermost ³P₀ level and thus the ${I}^{145}/{I}^{63}$ ratio. The relatively low line intensity ratio at high volume densities and low column densities is increased by the trapping resulting from increased column density. As the trapping increases, the excitation temperature of the 63 μm line begins to rise toward the kinetic temperature, while that of the 145 μm line first drops to quite low values (see Figure 4) but then also increases. The result is that, for a narrow range of column densities, the line ratio is almost independent of volume density. For higher column densities, the intensity ratio increases with increasing volume density up to a maximum value and then drops for higher n(H₂). This is seen most clearly for T_k = 400 K.

The ratio of the antenna temperatures, as discussed in Section 2.7.1 and described by Equation (31), is converging to a single value independent of n(H₂) and N(O^o). The behavior of the intensity ratio is not so easily calculated, due to the different line broadening due to different optical depths of the two transitions (see, e.g., Figure 7). The beginnings of this trend are seen especially clearly in the right-hand side of Figure 12 for T_k = 400 K.

4.2. Comparison with Other Models and Analysis of Observations

4.2.1. Molpop–CEP Output and Comparison with Results from LVG Calculation

The intensities of the [O I] fine structure lines alone can be presented in many different ways; we here give results in one format that is in effect complementary to that used previously. With only two observed quantities (two line intensities or one intensity and one intensity ratio), we cannot uniquely solve for the three key parameters: O⁰ column density, H₂ density, and kinetic temperature. Here, we fix the kinetic temperature and show the results in I(63 μm)–I(145 μm)/I(63 μm) space for different H₂ densities and O⁰ column densities. Figures 13–15 show the results for H₂ densities equal to 10³, 10⁴, and 10⁵ cm⁻³, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each figure, we plot the line intensity ratio versus the 63 μm line intensity at a fixed kinetic temperature for different atomic oxygen column densities. The points on each curve denote O⁰ column densities increasing in steps of 10^0.25 = 1.78, with the large squares indicating N = 1 × 10¹⁸ cm⁻².

The black curves show the results from the RADEX LVG code (van der Tak et al. 2007) for a kinetic temperature of 100 K. These curves are generally similar in form to those from Molpop–CEP, but have significantly higher line ratios for the same 63 μm line intensities at the same temperature. The RADEX calculation utilizes older deexcitation rate coefficients included in that program's database that result in higher line ratios, due to the relatively larger excitation rate to the upper, ³P₀, level of the triplet. The RADEX results thus mimic the results from Molpop–CEP at higher kinetic temperatures, conditions which also enhance the rate of collisions producing a ³P₂ → ³P₀ transition.

4.2.2. Comparison with PDR Models

The Molpop–CEP results can only be compared approximately with those from model PDR calculations such as those presented by Kaufman et al. (1999), inasmuch as the latter models a slab heated by an external UV field. Thus, while the density is uniform, the temperature drops substantially as we move inward from the surface. The form of the hydrogen is calculated based on a photochemical model, so it varies from an atomic form at the surface exposed to UV to largely molecular at extinctions ≥1 mag.

For a density of 10⁴ cm⁻³, Kaufman et al. (1999) determine a surface temperature of 100 K for G₀ = 10², log I(63 μm) = −3.4, and log I(145 μm)/I(63 μm) = −1.33. This is close to the value we obtain here for the same density, N(O⁰) = 1.8 × 10¹⁷ cm⁻², and T = 250 K. The ratio I(145 μm)/I(63 μm) increases with kinetic temperature for subthermal excitation, as discussed in Section 2.2.1 and illustrated in Figure 5. The gas temperature in the PDR model deceases from its surface value of 100 K, so the higher temperature required for the Molpop–CEP uniform-temperature slab compared to the PDR model cannot be a consequence of the temperature structure in the latter. The difference is likely a combination of the different collision rate coefficients and the change from molecular to atomic hydrogen as one moves toward the illuminated surface of the PDR.

4.2.3. Comparison with Observations of the Carina Nebula and Presence of Foreground Absorption

Figures 13–15 include intensities and intensity ratios from Oberst et al. (2011). The data presented in that paper included the [O I] observations of Mizutani et al. (2004), who observed both [O I] transitions in the central portion of the Carina Nebula. These observations did not resolve the lines spectrally, but having both transitions observed systematically offers a valuable point for evaluating the modeling tool described here. Mizutani et al. (2004) noted some positions with very high values of I(145 μm)/I(63 μm), extending up to 0.09. While details of their modeling are limited, they concluded that these values can only be obtained for very high temperatures, consistent with Kaufman et al. (1999) and with the results presented here. Our modeling does not extend to the required ≃1000 K temperatures, but Mizutani et al. (2004) claim these are excluded in any case by the ratio of [C ii] to [O I].

An alternative explanation of these high ratios is that they are due to a reduction in the 63 μm intensity that is due to absorption by lower-excitation foreground material (see, e.g., Saraceno et al. 1998). If we consider 50 K material having density 10⁴ cm⁻³ (either a low-density portion of the PDR or unrelated foreground material), the excitation temperature of the 63 μm line will be about 26 K (see Figure 2 for T_k = 100 K), and the fraction of atomic oxygen in the ³P₁ state ≃10⁻⁴. Consequently, the optical depth of the 145 μm line will be dramatically less than that of the 63 μm line. Using Equation (23) and values of τ₀ given in Table 3, for N(O⁰) = 6.4 × 10¹⁷ cm⁻² and Δv = 5 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, τ⁶³ = 0.65, while τ¹⁴⁵ ≃ 10⁻⁴. Consider as the "background source" a position with relatively strong PDR emission such as those observed by Mizutani et al. (2004) and Oberst et al. (2011), with I⁶³ = 1 × 10⁻⁴ erg s⁻¹ cm⁻² sr⁻¹ and I¹⁴⁵/I⁶³ = 0.044. With the above foreground optical depths for the two [O I] transitions, we see that the 63 μm intensity is reduced by a factor of ≃2 (including the very modest emission from the foreground cloud), while that of the 145 μm line is essentially unaffected. The result is that the intensity ratio is increased from 0.044 to 0.085.

Thus, a pair of line intensities that would be very difficult to fit now suggest PDR emission that is satisfactorily represented by a PDR with T_k ≃ 250 K, n(H₂) ≃ 10⁵ cm⁻³, and N(O⁰) ≃ 1 × 10¹⁶ cm⁻², with the less-excited, moderately optically thick foreground component discussed above. These parameters are not unique, but indicate that correcting for the foreground absorption can result in reasonable parameters for the atomic oxygen emission from the PDR instead of the observational results being problematic. The presence of foreground absorption can possibly be traced by the line profile, and this is another important aspect in which velocity-resolved spectra can be extremely valuable.

We have presented results from a step in developing more complete and realistic models with which to interpret observations of the submillimeter fine structure transitions of atomic oxygen. While observations of the [O I] 63 μm line, especially those that spectrally resolve the emission, have been limited, and those of the 145 μm line even less plentiful, observations of this atom have already proven to be powerful diagnostics of Galactic PDRs and widely employed as a tracer of star formation in Galactic sources and external galaxies. One significant development has been to go beyond the upper frequency limit of the HIFI instrument on Herschel (de Graauw et al. 2010), which now enables SOFIA observations of both [O I] lines by upGREAT (Risacher et al. 2016).

Observation of velocity-resolved spectra will give significantly more ability to interpret the complex structure of the PDR regions producing emission, which are often accompanied by less-excited [O I] associated with the emission region that produces self-absorption features, and also unrelated foreground [O I] responsible for additional absorption features. We have thus utilized the Molpop–CEP code to study [O I] emission. This radiative transfer and statistical equilibrium code accounts for photon trapping and produces emergent line profiles that can readily be compared with observations. We here present results from a uniform slab model that can be used as a first approximation in analyzing observations.

We present some salient features of the collisional excitation of the [O I] fine structure levels based on recently published quantum calculations by Lique et al. (2018). We confirm earlier analyses that found population inversion of the shorter-wavelength 145 μm transition, and we find that radiative trapping of the much more optically thick lower 63 μm transition suppresses this inversion. We find that the 63 μm line is often optically thick and that trapping increases its excitation temperature if the density is below the critical density, 5 × 10⁵ cm⁻³ and 8 × 10⁵ cm⁻³ for H₂ and H⁰, respectively. The less-excited [O I] at the surface of the slab readily produces self-absorption features in the 63 μm line profiles, similar to those often observed.

However, the absorption features produced by a cloud modeled as a slab with uniform physical conditions are not nearly as deep as are sometimes observed, suggesting that more complete models incorporating variations in temperature and density within the PDR are necessary. The temperature gradients are almost certainly present, due to heating by an external source, and the density variations are reasonable as well as being suggested by observation of different PDR tracers (e.g., Mookerjea et al. 2019).

We gratefully acknowledge major assistance from Moshe Elitzur and Andrés Asensio Ramos in using and understanding details of the Molpop–CEP code. We thank François Lique for very helpful comments about [O I] excitation and David Flower for explanation of collisional selection rules. We thank the anonymous reviewer for constructive comments that improved this paper. This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

For generality, we here label the three fine structure levels 1, 2, and 3 in order of increasing energy, denoting by n_i the volume density of the species of interest in level i. The individual rates of population transfer from level i to level j are written R_ij and include radiative and collisional rates. The rate equations for a steady state can be solved to yield two expressions for the ratios of populations, which are given as Equations (33) and (34):

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{3}}{{n}_{2}}=\displaystyle \frac{({R}_{21}+{R}_{23})({R}_{12}+{R}_{13})-{R}_{21}{R}_{12}}{{R}_{31}{R}_{12}+{R}_{32}({R}_{12}+{R}_{13})}\end{eqnarray} \tag{ 33 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{2}}{{n}_{1}}=\displaystyle \frac{{R}_{13}{R}_{32}+{R}_{12}({R}_{31}+{R}_{32)}}{({R}_{21}+{R}_{23})({R}_{31}+{R}_{32})-{R}_{23}{R}_{32}}.\end{eqnarray} \tag{ 34 }$$

These are convenient for analysis of the behavior of the line ratio in particular cases (e.g., Section 2.2). The constraint equation for the total density n_t of the species of interest is n_t = n₁ + n₂ + n₃, which we write to give n₁ using the above ratios as²

$$\begin{eqnarray}&&{n}_{1}={n}_{t}{\left[1+\displaystyle \frac{{n}_{2}}{{n}_{1}}+\displaystyle \frac{{n}_{3}}{{n}_{2}}\displaystyle \frac{{n}_{2}}{{n}_{1}}\right]}^{-1}.\end{eqnarray} \tag{ 35 }$$

The absolute values of the level populations are then simply obtained by multiplication of the preceding expressions.

While complete treatments include an external radiation field, and radiative trapping where appropriate, we can gain an appreciation of important aspects of behavior by considering that radiative rates are only the downward rates that are due to spontaneous emission. Collision rates are generally given as downward rate coefficients and multiplied by the density of colliding particles to get the rates. The upward collision rates are determined from a detailed balance.