High-resolution SOFIA/EXES Spectroscopy of Water Absorption Lines in the Massive Young Binary W3 IRS 5

Massive stars reach the main sequence while deeply embedded and still accreting. While rich chemistry is driven as the central object warms and ionizes the environment, complicated physical activities such as accretion disks, outflows, shocks, and disk winds are involved (Beuther et al. 2007; Cesaroni et al. 2007; Zinnecker & Yorke 2007). However, the large distances to the observers, the high extinction at optical and near-infrared wavelengths, and the highly clustered environment impede a clear understanding of their formation and evolution processes.

High spectral resolution, pencil beam absorption line studies at mid-infrared (MIR) wavelengths provide a unique opportunity to probe the embedded phases in massive star formation (Lacy 2013). First, the MIR continuum originates from the photosphere of the disk at a distance of tens to a few hundreds of astronomical unit (Beltrán & de Wit 2016; Frost et al. 2021). The effective pencil beam, therefore, avoids beam dilution issues that submillimeter observations are subject to. Second, molecules without dipole moments such as C₂H₂ and CH₄, which are among the most abundant carbon-bearing molecules, can only be observed through their rovibrational spectra in the infrared. Therefore, MIR spectroscopy at high resolution is a critical tracer of the physical conditions, the chemical inventory, and the kinematics of structures close to the massive protostars.

Among the rich chemical inventory in the regions associated with the protostars, water is of fundamental importance because it is one of the most abundant molecules in both the gas and ice phases. As its abundance varies largely between warm and cold gas (e.g., Draine et al. 1983; Kaufman & Neufeld 1996; Bergin et al. 2002), water has a powerful diagnostic capability in probing physical conditions (van Dishoeck et al. 2021). However, due to its prevalence in the Earth's atmosphere, water is very difficult to observe from the ground. We present in this work the power of combining the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012) that observes between 39,000 and 45,000 ft and the Echelon Cross Echelle Spectrograph (EXES; Richter et al. 2018) spectrometer, which resolves lines to several km s⁻¹ to make the most of the diagnostic capability of water.

Previous studies of the MIR water absorption spectrum toward massive protostars (Boonman & van Dishoeck 2003; Boonman et al. 2003) have revealed a rich spectral content. High spectral resolution observations can provide much insight into the characteristics of these regions (e.g., Knez et al. 2009; Indriolo et al. 2015; Barr et al. 2018; Dungee et al. 2018; Indriolo et al. 2018; Rangwala et al. 2018; Goto et al.2019; Barr et al. 2020; Indriolo et al. 2020; Nickerson et al. 2021; Barr et al. 2022a, 2022b; Li et al. 2022; Nickerson et al. 2023). Water has also been studied at high spectral resolution via pure rotational transitions, using the heterodyne instruments on board SWAS, Odin, and Herschel Space Observatory (e.g., Snell et al. 2000; Wilson et al. 2003; Chavarría et al. 2010; Karska et al. 2014). Because of their limited spatial resolution, these observations mainly probed the large-scale environment of these sources.

We conducted high spectral resolution (R = 50,000) spectroscopy from 5–8 μm with EXES on board SOFIA toward the hot core region close to the massive binary protostar W3 IRS 5, which is a luminous, massive protostellar binary in transition from the embedded to the exposed phase of star formation (van der Tak et al. 2000, 2005). W3 IRS 5 is oriented along the northeast-to-southwest direction with a separation of 12 (∼2800 au at ${2.3}_{-0.16}^{+0.19}$ kpc; Navarete et al. 2019). Following the nomenclature in van der Tak et al. (2005), we refer to the northeastern object of the binary as MIR1 and the southwestern one as MIR2.

Millimeter studies reveal complex structures surrounding the binary on a large scale (10³–10⁵ au), including a hot molecular core with outflows, jet lobes, shocks, and a circumbinary toroid (Imai et al. 2000; Rodón et al. 2008; Wang et al. 2012, 2013; Purser et al. 2021). In the infrared M-band (4.7 μm), high spectral resolution (R = 88,100) observations of rovibrational transitions of CO and its isotopologues (Li et al. 2022) show various absorbing structures along the pencil beam line of weight against the pencil beam MIR background. Physical conditions such as temperature and column density have been derived from these CO observations. While MIR1 and MIR2 are spatially resolved in Li et al. (2022), the absorbing components are identified as a shared envelope (∼50 K at −38.5 km s⁻¹), several foreground clumps produced by either J- or C-shocks (200–300 K from −100 to −60 km s⁻¹), and blobs that are likely associated with the circumstellar disks (∼500 K from −55 to −38.5 km s⁻¹). The M-band study sets up the backdrop for understanding the origin of spatially unresolved water absorption lines in this paper, and we refer readers to Table 7 in Li et al. (2022) for a complete list of properties of each decomposed CO component.

This paper presents the rich spectrum of rovibrational (the ν₂ band) water lines observed in absorption toward W3 IRS 5. The high spectral resolution spectroscopy allows us to separate and identify individual velocity components that are linked to different stars in the W3 IRS 5 binary and to derive the temperatures, the level-specific column densities, as well as the total column densities (and/or the abundances). We describe our observations and data reduction in Section 2 and our analysis methods with both the rotation diagram and the curve-of-growth in Section 3. We present in Section 4 the properties of multiple dynamical components and set up the connection between these water components to CO components that are identified in Li et al. (2022). We discuss the chemical abundances along the line of sight based on the CO-to-H₂O connection in Section 5.

Object W3 IRS 5 was observed with the EXES spectrometer (Richter et al. 2018) on board the SOFIA observatory from 2021 June to 2022 February as part of programs 08_0136 and 09_0072. Archival data observed in programs 07_0063 and 76_0004 were also included in this study. The full spectral survey covers a wavelength range of 5.3–7.9 μm in 18 observational settings and was observed under the HIGH-LOW cross-dispersed mode. The observational parameters of all 18 settings are listed in Table 1. For all settings, the slit width is fixed to 32 to limit slit losses perpendicular to the slit at an SOFIA point-spread function (PSF) FWHM (see ∼30–35), which provides a spectral resolution of R ∼ 50,000, or equivalently, a velocity resolution of 6 km s⁻¹. The slit length is dependent on the wavelength and the angle of the echelle grating and is in the range of 136–343 after accounting for the anamorphic magnification (Figure 1). Off-slit nodding was applied to remove the background night sky emission and the telescope thermal emission.

Zoom In Zoom Out Reset image size

The EXES data were reduced with the SOFIA Redux pipeline (Clarke et al. 2015), which has incorporated routines originally developed for the Texas Echelon Cross Echelle Spectrograph (Lacy et al. 2003). The science frames were de-spiked and sequential nod positions subtracted, to remove telluric emission lines and telescope/system thermal emission. An internal blackbody source was observed for flat-fielding and flux calibration and then the data were rectified, aligning the spatial and spectral dimensions. The wavenumber solution was calibrated using sky emission spectra produced for each setting by omitting the nod-subtraction step. We used wavenumber values from HITRAN (Rothman et al. 2013) to set the wavelength scale. The resulting wavelength solutions have 1σ errors of 0.5 km s⁻¹.

Before measuring intrinsic lines in the source, the spectra had to be corrected for telluric absorption. The telluric absorption lines can be separated from lines that are from W3 IRS 5 thanks to the Doppler shifts (Table 1). Ideally, removing telluric absorption is done by taking a spectrum of a featureless hot star immediately before or after the target observation. This method also has the advantage of removing instrumental baseline effects, such as fringing, present in both spectra. Unfortunately, it was impractical to schedule calibrators for every part of the W3 IRS 5 survey due to the limitations of airborne observation scheduling and the scarcity of bright standard stars. Only one survey setting was observed with the adjacent calibrator star Sirius (7.28–7.46 μm; see Table 1). Thus, we relied on atmospheric transmittance models created with the Planetary Spectrum Generator (PSG; Villanueva et al. 2018) and tuned the H₂O column and pressure for the longitude, latitude, and altitude of the observations (Appendix A) to remove telluric features, including minor features that may overlap with water lines from W3 IRS 5. The quality of these corrections was verified by telluric lines not overlapping with W3IRS5 lines. To establish the reliability of the construction procedure of the PSG models, we compared the Sirius spectrum with the PSG model built from 7.28–7.46 μm (see Appendix B). We conclude that the telluric correction yields uncertainties on the equivalent widths at the level of 10%. After the PSG models are divided, further data reduction such as removing the local baseline is still conducted.

Since the distance between the binary (12) is smaller than the SOFIA PSF (∼3''–35), W3 IRS 5 is not spatially resolved in the observed spectra. As both sources contribute to the observed flux in our slit, the derived equivalent width of an absorption associated with one IR source will be diminished by the continuum emission from the other source. The two sources have very comparable MIR brightness (see the SpeX observations in Li et al. 2022), and we include a factor of 2 in our analysis to account for this effect throughout this study unless otherwise stated. In other words, we think the relative absorption depth should be two times deeper because the continuum level should be two times weaker. It should be mentioned that the actual contribution of the nonabsorbing source will depend on the slit orientation over the source, and this changes slightly between the different grating settings (see Figure 1). We have decided to accept this additional source of uncertainty to the results in view of the uncertain brightness distribution of the two sources on the sky.

By comparing with the spectra constructed via the existing laboratory line information (HITRAN; Rothman et al. 2013) and LTE models, we identified over 180 H₂O ν₂ = 1–0 and about 90 H₂O ν₂ = 2–1 absorption lines in this survey. As shown in Figure 2, the velocity ranges of the absorbing components are very similar to those present in ¹³CO or in high-J ¹²CO lines (Li et al. 2022). All components are blueshifted compared to or are located at the cloud velocity v_LSR= −38 km s⁻¹ (van der Tak et al. 2000). Each kinematic component is characterized by different excitation conditions.

Zoom In Zoom Out Reset image size

We fit the ν₂ = 1–0 absorption line profiles with a sum of multiple Gaussians indicated as the ith component in the form of

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\mathrm{obs}}(v)}{{I}_{c}}=1-\displaystyle \sum _{i}\displaystyle \frac{{W}_{i}}{{\sigma }_{v,i}\sqrt{2\pi }}\exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\left(v-{v}_{\mathrm{LSR},i}\right)}^{2}}{{\sigma }_{v,i}^{2}}\right)\end{eqnarray} \tag{ 1 }$$

using the curve_fit function in scipy. In Equation (1), I_obs(v) is the observed line intensity profile in the velocity space, I_c is the intensity of the continuum, W_i is the area or the equivalent depth of the Gaussian component, σ_v,i is the standard deviation, and v_cen,i is the center velocity of the component. The fitting results are determined by restricting the range of the free parameters. Specifically, two components centered at ∼ − 39.5 km s⁻¹ and −54.5 km s⁻¹ are present across all energy levels, while the component centered at −45 km s⁻¹ only appears in transitions originating from relatively low-energy states (<800 K). As we will explain in Section 3.1, we name the three components at ∼ − 39.5 km s⁻¹, −45 km s⁻¹, and −54.5 km s⁻¹ as "H1," "W," and "H2." For very-low-energy states (<200 K), an extra component at −39.5 km s⁻¹ is possibly needed for a better fitting result (see Appendix C), and we discuss its existence and the implications in Section 5.1.

We therefore fit two and three Gaussians above and below 800 K, separately. Above 800 K, we constrain the velocity center to −43 to −36 km s⁻¹ for the "H1" component, and the "H2" component with a relative velocity of −17 to −14 km s⁻¹ guided by an initial guess. The line widths, σ_v , were constrained to 4–6 and 7–9 km s⁻¹. We apply the same constraints to the ν₂ = 2–1 absorption line profiles. On the other hand, we think the profiles of the low-energy transitions (<800 K) are a blend of more than two components, including the rather weak "W" and "H2" components. After a careful assessment, we decided to limit the number of free parameters in the fit to these transitions as, from a physical perspective, we expect that transitions from a given temperature component will occur at consistent central velocities. We keep using a range for the velocity dispersion in the fitting procedure as we did for the high-energy transitions. More specifically, for the low-energy transitions, we fix the rather weak "W" and "H2" with the central velocity, v_LSR, of −45 and −54.5 km s⁻¹, and the velocity dispersion parameter, σ_v , from 4–5 km s⁻¹. We note that the H2 (−54.5 km s⁻¹) component disappears in some of the low-energy transitions (see the second panel in Figure 2) as a result of the opacity effect (see Section 3.1). In this case, we only fit the line profile with one Gaussian for the "H1" component. We list in Appendix E in detail the parameters of each individual component such as the central velocity, velocity width, and equivalent widths, and report "H1" as −39.1 ± 2.4 km s⁻¹, and "H2" (>800 K) as −54.1 ± 3.0 km s⁻¹, respectively.

3.1. Rotation Diagram Analysis and the Opacity Problem

After the distinct velocity components are determined, we use a rotation diagram to provide a first view of their characteristics. If the absorption lines are optically thin, we can get the column density N_l in the lower state of a transition directly from the integrated line profile by

$$\begin{eqnarray}&&{N}_{l}=8\pi /({A}_{{ul}}{\lambda }^{3}){g}_{l}/{g}_{u}\int \tau (v){dv},\end{eqnarray} \tag{ 2 }$$

in which λ is the wavelength, A_ul is the Einstein A coefficient, g_l and g_u are the statistical weight of the lower and upper levels, v is the velocity, and

$$\begin{eqnarray}&&\tau (v)=-\mathrm{ln}({I}_{v}/{I}_{c}),\end{eqnarray} \tag{ 3 }$$

where I_v and I_c are the intensity of the absorption line and the continuum, and I_v of each identified velocity component in a given transition is derived from the Gaussian fitting parameters (Equation (1)). All spectral line parameters used in this study are adopted from HITRAN (Rothman et al. 2013).

If the foreground absorbing gas is in LTE, the population of one rotational level can be described with the Boltzmann equation (Goldsmith & Langer 1999),

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{l}}{{g}_{l}}=\displaystyle \frac{{N}_{\mathrm{tot}}}{Q({T}_{\mathrm{ex}})}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{E}_{l}}{{k}_{B}{T}_{\mathrm{ex}}}\right),\end{eqnarray} \tag{ 4 }$$

in which T_ex is the excitation temperature, N_tot is the total column density, ΔE_l is the relative energy between the excitation state and the ground state energy of the vibration state, and is equal to E_l for ν₂ = 1–0 transitions. ¹² Q(T_ex) is the partition function, and k_B is the Boltzmann constant. Specifically, ln(N_l /g_l ) and E_l /k_B of all absorption lines constructs the so-called rotation diagram. For a uniform excitation temperature, ln(N_l /g_l ) and E_l /k_B fall on a straight line. The inverse of the slope represents the temperature, and the intercept represents the total column density over the partition function.

We present the rotation diagrams of the rovibrational H₂O lines from W3 IRS 5 in Figure 3. For the three velocity components centered at −39.5, − 45, and −54.5 km s⁻¹ identified in ν₂ = 1–0 transitions, the derived temperatures are 807 ± 58, 200 ± 18, and 669 ± 57 K (Table 2). The three components are hence named "H1," "W," and "H2," in which "W" and "H" stand for "warm" and "hot," respectively, following the nomenclature in Li et al. (2022). For ν₂ = 2–1 transitions, only the components "H1" and "H2" at −39.5 and −54.5 km s⁻¹ are identified, and the derived temperatures are 703 ± 60 and 946 ± 170 K (Table 2). ¹³

Zoom In Zoom Out Reset image size

Table 2. Physical Conditions of Decomposed Components

Component	"W," −45 km s−1	"H1," −39.5 km s−1		"H2," −54.5 km s−1
Transitions	ν2 = 1–0	ν2 = 1–0	ν2 = 2–1	ν2 = 1–0	ν2 = 2–1
Rotation Diagram
Tex (K)	200 ± 18	807 ± 58	703 ± 60	669 ± 57	946 ± 170
Ntot ( cm−2)	(6.6 ± 4.0) × 1017 a	(1.2 ± 0.6) × 1018	(7.9 ± 3.7) × 1016	(6.2 ± 3.4) × 1017	(8.1 ± 5.6) × 1016
Slab Model
fc	⋯	0.4	0.4	0.3	⋯
b ( km s−1)	⋯	2.8	2.8	3.5	⋯
Tex (K)	⋯	${471}_{-15}^{+14}$	${654}_{-191}^{+135}$	${600}_{-27}^{+28}$	⋯
Ntot ( cm−2)	⋯	${2.5}_{-0.2}^{+0.3}\times {10}^{19}$	${2.8}_{-0.7}^{+0.5}\times {10}^{17}$	${5.3}_{-0.6}^{+0.3}\times {10}^{18}$	⋯
${\chi }_{r,{\min }}^{2}$	..	1.7	0.65	3.8	⋯
Disk Model b
b ( km s−1)	⋯	2.8	2.8	3.5	⋯
Tex	⋯	${491}_{-14}^{+13}$	${691}_{-212}^{+122}$	${612}_{-30}^{+27}$	⋯
Abun. (w.r.t H)	⋯	${2.6}_{-0.2}^{+0.1}\times {10}^{-3}$	${5.0}_{-3.5}^{+1.4}\times {10}^{-5}$	${5.1}_{-0.5}^{+0.4}\times {10}^{-4}$	⋯
${\chi }_{r,{\min }}^{2}$	⋯	1.7	0.34	2.3	⋯

Notes.

^a This value is corrected by assuming that "W" covers protobinary and that the absorption intensity is not diluted (see Section 4.2). ^b Pure absorption with parameter of 1 is assumed in this table. See Table 3 for results with an of 0 and 0.5.

Download table as:  ASCIITypeset image

The derived total column densities of "H1," "W," and "H2" in the ν₂ = 0 state are (1.2 ± 0.6) × 10¹⁸, (3.3 ± 2.0) × 10¹⁷, and (6.2 ± 3.4) × 10¹⁷ cm⁻², respectively (Table 2). However, it is noticeable that the large scatter in the rotation diagram exceeds what the error bars can account for in components H1 and H2. Specifically, as shown in Figure 3, for transitions that share the same lower state (and the same g_l ), the difference in the derived N_l is up to an order of magnitude.

We can further illustrate this problem by directly comparing the line profiles of the transitions that have a common lower energy level, and present such a comparison in Figures 2 and 4. As shown in Figures 2 and 4, lines that share the same lower level but different oscillator strength (f_ul ) or Einstein A do not necessarily have the same absorption intensity and/or the equivalent width, which is defined as

$$\begin{eqnarray}&&{W}_{\nu }=\int (1-{I}_{\nu }/{I}_{c})d\nu \end{eqnarray} \tag{ 5 }$$

in the frequency space. This behavior where transitions with high values of Einstein A fall systematically below the relation provided by transitions with low values of Einstein A in the rotation diagram (Figure 3) is characteristic for opacity effects. Specifically, when transitions are optically thick, an increase in absorption strength (due to an increase in Einstein A) results in only a small increase in line width (not in depth) and hence marginally increases the equivalent width. This effect was noted earlier in Indriolo et al. (2020) and Barr et al. (2022a) in studies of MIR H₂O rovibrational lines toward massive protostars AFGL 2136 and AFGL 2591.

Zoom In Zoom Out Reset image size

3.2. Curve-of-growth Analysis

Considering that τ(v) in Equation (2) does not represent an optically thin Gaussian core, the definition of the equivalent width in Equation (5) can be written as

$$\begin{eqnarray}\begin{array}{rcl}{W}_{\nu } & = & \displaystyle \int (1-{I}_{\nu }/{I}_{c})d\nu \\ & = & \displaystyle \int (1-{{\rm{e}}}^{-\tau (\nu )})d\nu \\ & = & \displaystyle \int (1-{{\rm{e}}}^{-{\tau }_{p}H(a,v)})d\nu \end{array}\end{eqnarray} \tag{ 6 }$$

in which the integration is over the frequency ν. The Voigt profile H(a, v) is defined as (Equations (9)–(34) in Mihalas1978):

$$\begin{eqnarray}&&H(a,v)=\displaystyle \frac{a}{\pi }{\int }_{-\infty }^{+\infty }\displaystyle \frac{{e}^{-{y}^{2}}{dy}}{{\left(v-y\right)}^{2}+{a}^{2}}.\end{eqnarray} \tag{ 7 }$$

The parameter, v, is defined as

$$\begin{eqnarray}&&v=\displaystyle \frac{\nu -{\nu }_{0}}{{\rm{\Delta }}{\nu }_{D}}\end{eqnarray} \tag{ 8 }$$

and represents the shift from the line center in Doppler units. Δν_D is the Doppler width in frequency space. The parameter a is the damping factor. The peak optical depth τ_p is

$$\begin{eqnarray}&&{\tau }_{p}=\displaystyle \frac{\sqrt{\pi }{e}^{2}}{{m}_{e}{bc}}{N}_{l}{f}_{{lu}}\lambda .\end{eqnarray} \tag{ 9 }$$

In Equation (9), e is the electron charge, m_e is the electron mass, c is the speed of light, and f_lu is the oscillator strength. The Doppler parameter in velocity space, b, is related to the FWHM of an optically thin line by ${\rm{\Delta }}{v}_{\mathrm{FWHM}}=2\sqrt{\mathrm{ln}2}b$.

Equation (9) clearly shows that for lines that share the same lower level and have the same N_l , a difference in f_lu λ will lead to different τ_p . Defining log ₁₀(f_lu λ) as the representative for the line strength, lines with a larger line strength have larger equivalent widths (as shown in Figure 4), and as a result, the N_l derived via Equation (2) will be underestimated (as shown in Figure 3).

3.2.1. Slab Model of a Foreground Cloud

A curve-of-growth analysis (Rodgers & Williams 1974) is required to reconcile the opacity problem and to correctly derive N_l and N_tot (Equation (9.27) in Draine 2011):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{W}_{\lambda }}{\lambda }\approx \\ \left\{\begin{array}{ll}\displaystyle \frac{\sqrt{\pi }b}{c}\displaystyle \frac{{\tau }_{p}}{1+{\tau }_{p}/(2\sqrt{2})} & \mathrm{for}\,\tau \lt 1.254\\ \displaystyle \frac{2b}{c}\sqrt{\mathrm{ln}({\tau }_{p}/\mathrm{ln}2)+\displaystyle \frac{\gamma \lambda }{4b\sqrt{\pi }}({\tau }_{p}-1.254)} & \mathrm{for}\,\tau \gt 1.254\end{array}\right..\end{array}\end{eqnarray} \tag{ 10 }$$

We specifically note that this correction applies to an absorbing foreground slab model. In the equations above, the definitions of all parameters follow Equation (9). The parameter γ is the damping constant of the Lorentzian profile and is of the order of 10 for radiative damping. We stress that for H₂O lines discussed in this paper, the Lorentzian line width that γ corresponds to is 10⁻⁹ km s⁻¹, and is negligible compared to the observed Doppler width (a few kilometeres per second).

For such an absorbing slab model, the emission from the foreground is negligible against the representative background temperature of ∼600 K in this study (see Li et al. 2022). Furthermore, if the foreground cloud does not cover the entire observing beam, a covering factor f_c (0 ≤ f_c ≤ 1) has to be accounted for, and Equation (3) is modified to:

$$\begin{eqnarray}&&{I}_{v}={I}_{{\rm{c}}}(1-{f}_{c}(1-{{\rm{e}}}^{-\tau (v)})),\end{eqnarray} \tag{ 11 }$$

and the left-hand side of the Equation (10) is modified to W_λ /(λ f_c ).

3.2.2. Stellar Atmosphere Model of a Circumstellar Disk

The absorption can also occur in an accretion disk scenario in the system of a forming massive star (Barr et al. 2020, 2022a; Li et al. 2022). Specifically, the disk has an outward-decreasing temperature gradient from the mid-plane to the surface. Such disks show absorption lines because the thermal continuum from the dust is mixed with the molecular gas. For such a scenario, we adopt the curve-of-growth of the stellar atmosphere model in which the continuum and the line opacities are coupled. The residual flux,

$$\begin{eqnarray}&&{R}_{\nu }\equiv {I}_{\nu }/{I}_{c},\end{eqnarray} \tag{ 12 }$$

where I_ν is the intensity of the absorption line, can then be approximated by the Milne-Eddington model (Mihalas 1978) which assumes a gray atmosphere. The absorption line profile may originate in pure absorption or scattering. The parameter characterizes the line formation, and can take the form of 1 (pure absorption), 0 (pure scattering), or between 0 and 1 (a combination of scattering and absorption). We refer to Appendix A in Barr et al. (2020) for details of the line residual flux expected from this model.

The curve of growth in the stellar atmosphere model is constructed via the equivalent width versus β₀, the ratio of the line opacity at the line center, κ_L (ν = ν₀), to the continuum opacity κ_c :

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{W}_{\nu }}{2{\rm{\Delta }}{\nu }_{D}} & = & {\displaystyle \int }_{0}^{+\infty }(1-{R}_{v}){dv}\\ & = & {A}_{0}{\displaystyle \int }_{0}^{+\infty }{\beta }_{0}H(a,v){[1+{\beta }_{0}H(a,v)]}^{-1}{dv},\end{array}\end{eqnarray} \tag{ 13 }$$

in which

$$\begin{eqnarray}\begin{array}{rcl}{\beta }_{0} & = & \displaystyle \frac{{\kappa }_{L}(\nu ={\nu }_{0})}{{\kappa }_{c}}\\ & = & \displaystyle \frac{{A}_{{ul}}{\lambda }^{3}}{8\pi \sqrt{2\pi }{\sigma }_{v}}\displaystyle \frac{{g}_{u}}{{g}_{l}}\displaystyle \frac{{N}_{l}}{{\sigma }_{c}{N}_{H}}\left(1-\displaystyle \frac{{g}_{l}}{{g}_{u}}\displaystyle \frac{{N}_{u}}{{N}_{l}}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

where we have used the continuum opacity, κ_c , equal σ_c N_H where σ_c is the dust cross section per H atom, H(a, v) is the Voigt function that gives the line profile in velocity space, and v is defined as the velocity shift with respect to the line center in units of the Doppler width (see Equation (7)). The damping factor a = γ λ/b is of the order of 10⁻⁸ for H₂O rovibrational lines. The parameter A₀ is the central depth of an opaque line. Its exact value is determined by the radiative transfer model of the surface of the disk and is related to the gradient of the Planck function. For a gray atmosphere and lines in pure absorption, A₀ is ∼0.5–0.9 from 900–100 K (see Appendix A in Barr et al. 2020). The dispersion in velocity space, σ_v , is transformed from the Doppler parameter, $b/\sqrt{2}$, with which we convert W_ν to the velocity space:

$$\begin{eqnarray}&&\displaystyle \frac{{W}_{\nu }}{2{\rm{\Delta }}{\nu }_{D}}=\displaystyle \frac{{W}_{v}}{2b}.\end{eqnarray} \tag{ 15 }$$

We can disregard the bracketed item in Equation (14) if stimulated emission is negligible.

We adopt a value of 7 × 10⁻²³ cm²/H-nucleus for σ_c following Barr et al. (2020), as it is appropriate for coagulated interstellar dust (Ormel et al. 2011). Theoretical fits to these curves of growth will provide abundances relative to the dust opacity. As this value of σ_c is adopted in the CO analysis (Li et al. 2022) as well, we are able to derive an absolute water-to-CO ratio once the CO correspondent component to water is identified via the kinematic information (see Section 4.2).

As shown in Section 3, the H₂O ν₂ = 1–0 lines are decomposed into three velocity components at −45, − 39.5, and −54.5 km s⁻¹, while the ν₂ = 2–1 absorption lines are decomposed into two components at −39.5 and −54.5 km s⁻¹. The components at the three different velocities, as explained in Section 3.1, are labeled as "W," "H1," and "H2" based on the preliminary estimation of their temperatures derived from the rotation diagram analysis under the optically thin assumption (Table 2).

We hereafter apply the curve-of-growth analyses to the ν₂ = 1–0 transitions from "H1 and 'H2" and ν₂ = 2–1 from "H1." The rotation diagrams of the three sets of transitions illustrate a large scatter indicative of the opacity problem discussed in Section 4.1, unlike for the other components and transitions. In Section 4.2, all of the water components are compared to and connected with warm and/or hot CO components. The implications of the vibrationally excited water lines are presented in Section 5.1.2.

4.1. Two Hot Physical Components: H1 and H2

We conduct the grid-search method (Li et al. 2022) on the (T_ex, N_tot) and (T_ex, abundance) grid in the curve-of-growth analyses for the slab model and the disk model, respectively. For the slab model, T_ex together with N_tot determines τ_p , or Nf_ul λ (Equation 10). We, therefore, search for the parameter combination of (T_ex, N_tot) that gives the smallest reduced χ² between the observed equivalent width W_λ (in the wavelength space) and that derived from the theoretical curve of growth (Equation (10)). Similarly, for the disk model, T_ex together with the abundance (N/N_H) determines β₀, the ratio of the line opacity to the continuum opacity. We search for the best (T_ex, abundance) combination that gives the smallest reduced χ² between the observed equivalent width W_ν (in the frequency space) and that derived from the theoretical curve of growth (Equation (13)).

When fitting the observed W_λ /λ to the theoretical curve-of-growth of a slab model, we note that the partial coverage, f_c , and the Doppler width, b (= $\sqrt{2}{\sigma }_{v}$), are degenerated parameters. Different combinations of f_c and b can provide as good fitting results, and we illustrate this point below for a similar case in the disk model. Observational results can put constraints on the range of f_c and b to some extent. For example, f_c are constrained as ∼0.4, which is twice the lower limit of the line intensities, ∼0.2. The factor of 2 is because of the dilution effect—the effect that the absorption line is against the continuum contributed by both MIR sources in W3 IRS 5 (see Section 2). As for the Doppler width b, its lower limit can be constrained by the thermal line width, and for gas of 500 K, ${\sigma }_{\mathrm{th}}=\sqrt{({k}_{{\rm{B}}}T)/(\mu {m}_{{\rm{H}}})}$ = 9.12 × 10⁻² km s⁻¹ (T/K)^0.5 μ^−0.5 ≈0.5 km s⁻¹. The upper limit of b can be constrained by the observed line width, which is b convolved with the instrument resolution ${\sigma }_{\mathrm{res}}$ = $c/(2\sqrt{2\mathrm{ln}2}R)$ = 2.5 km s⁻¹.

For the disk model, there is also a dependence on the chosen line width b, and the degree of absorption relative to the scattering in the line . We illustrate this in a more quantitative way. Take the "H1" component as an example, as presented in Table 3, while the results of the temperature and abundance depend on the adopted σ_v and , combinations of different σ_v and may provide comparable ${\chi }_{r,{\min }}^{2}$. Although we may have some control over σ_v , the value of is unconstrained. We, therefore, provide Table 3 as a reference for conditions when the lines are not due to pure absorption.

Table 3. Results of the Physical Component "H1" Derived from Different and σ_v in the Disk Model

	= 0			= 0.5			= 1
	${\chi }_{r,{\min }}^{2}$	Tex (K)	${X}_{{{\rm{H}}}_{2}{\rm{O}}}/{X}_{\mathrm{CO}}$	${\chi }_{r,{\min }}^{2}$	Tex (K)	${X}_{{{\rm{H}}}_{2}{\rm{O}}}/{X}_{\mathrm{CO}}$	${\chi }_{r,{\min }}^{2}$	Tex (K)	${X}_{{{\rm{H}}}_{2}{\rm{O}}}/{X}_{\mathrm{CO}}$
σv = 1.5 km s−1	2.02	579	1.19	1.79	449	2.45	1.88	419	3.71
σv = 2.0 km s−1	2.54	565	0.78	1.84	510	1.30	1.73	491	1.53
σv = 2.5 km s−1	3.14	598	0.61	1.99	560	0.90	1.88	521	1.20
σv = 3.0 km s−1	3.72	622	0.52	2.37	594	0.72	1.98	564	0.92

Note. (1) The 1σ uncertainty for the derived temperatures is of the order of ±10 K. (2) ${X}_{{{\rm{H}}}_{2}{\rm{O}}}/{X}_{\mathrm{CO}}$ is derived by assuming that "H1" and "MIR2-H1" in CO coexist.

Download table as:  ASCIITypeset image

For "H1" and "H2," we therefore choose the partial coverage, f_c , to 0.4 and 0.3, and set the Doppler width b (= $\sqrt{2}{\sigma }_{v}$) to 2.8 and 3.5 km s⁻¹, separately, since the chosen values provides the smallest ${\chi }_{r}^{2}$. For the ν₂ = 1–0 transition, we present the best-fitting results of the "H1" component for both the foreground and disk model in Figure 5 (results of "H2" are presented in Figure 12 in Appendix D) and summarize the derived properties in Table 2, assuming that the lines are formed in pure absorption ( = 0). In either the slab or the disk model, about half of the data points are located on the logarithmic part of the curve-of-growth, confirming that the corresponding absorption lines are optically thick. In the slab model, the curve-of-growth analysis "corrects" the underestimated total column densities in the rotation diagram by factors of 21 and 9 for "H1" and "H2," respectively. The derived temperatures are lowered to 471 K and 600 K, respectively (Table 2). In the disk model, as relative abundances are derived, the correction in column densities is quantified once the connection between CO and water components is established (see Section 4.2). The derived temperatures from the disk model are comparable to those derived from the slab model (Table 2). As for the ν₂ = 2–1 transition on "H1," the correction after the curve-of-growth analysis on the temperature and column density is insignificant (Figure 13 in Appendix D and Table 2), indicating that the optical depth effect is not severe.

Zoom In Zoom Out Reset image size

4.2. Connecting the Absorbing Components of H₂O and CO

High spectral resolution M-band absorption spectroscopy (R = 88,100) toward W3 IRS 5 at 4.7 μm revealed multiple kinematic components in ¹²CO and its isotopologues (Li et al. 2022). Specifically, IRTF/iSHELL spatially resolved the two protostars in W3 IRS 5. As CO lines also trace ambient gaseous components that absorb against the MIR background, building a connection between the components identified in H₂O and CO will help us to understand the origins of the gaseous H₂O components.

The velocity resolution of EXES (6 km s⁻¹) and iSHELL (3.4 km s⁻¹) enables the comparison between H₂O and CO components through their kinematic information. We note that we elected to link the observed H₂O and CO components through their similarity in velocity and temperature. The velocity width is not used for aiding the comparison because it is dependent on the optical depth and is degenerated with the column density (Table 3). We, therefore, present in Figure 6 the comparison between the absorption profiles and list in Table 4 the velocity and derived excitation temperature of components identified in the two species.

Zoom In Zoom Out Reset image size

Table 4. Comparing Water and CO Components

Comp.	Species	vLSR	Tex
		( km s−1)	(K)
"W"	H2O	−46	200 ± 18
MIR1-W1	CO	−43	${180}_{-14}^{+11}$
MIR2-W2	CO	−45.5	116 ± 7
"H1"	H2O (Slab)	−39.4 ± 2.1	${471}_{-15}^{+14}$
	H2O (Disk)	−39.4 ± 2.1	${491}_{-14}^{+13}$
MIR2-H1	CO	−37.5	${662}_{-28}^{+23}$
"H2"	H2O (Slab)	−54.1 ± 3.0	${600}_{-27}^{+28}$
	H2O (Disk)	−54.1 ± 3.0	${612}_{-30}^{+27}$
MIR2-H2	CO	−52	${482}_{-70}^{+97}$
MIR1-W1${}^{{\prime} }$	CO	−40, ∼ − 54	${709}_{-101}^{+136}$

Note. Physical conditions of CO components are adopted from Table 7 in Li et al. (2022). The adopted temperatures of hot CO components (MIR2-H1, MIR2-H2, MIR1-W1') are corrected values assuming a disk model.

Download table as:  ASCIITypeset image

As the different velocity components are characterized by different excitation temperatures, both the profiles of the H₂O and CO spectra are averaged in low-energy (<400 K) and high-energy (>400 K) transitions. As we have described in Section 3, we identified three components centered at −39.5, − 45, and −54.5 km s⁻¹ ("H1," "W," and "H2") in low-energy transitions, and two centered at −39.5 and −54.5 km s⁻¹ in high-energy transitions. We note that half of the transitions (6 of 12) with a −45 km/s component are at energies between 400 and 800 K.

In low-energy levels, both MIR1 and MIR2 have two components revealed in CO, with one centered at the v_sys = − 38 km s⁻¹ (MIR1-C1 and MIR2-C1 ¹⁴ ), and the other centered at −46 km s⁻¹ (MIR1-W1 and MIR2-W2). No low-energy CO components were found at −55 km s⁻¹. MIR1-C1 and MIR2-C1, as well as MIR1-W1 and MIR2-W2, have similar temperatures and column densities (see Table 7 in Li et al. 2022; so do temperatures of CO components mentioned below), and are regarded to cover MIR1 and MIR2 simultaneously.

The temperatures of MIR1-C1 and MIR2-C1 are low (∼50 K). MIR1-C1 and MIR2-C1 are therefore not related to any water components. In contrast, the temperatures of MIR1-W1 and MIR2-W2 are ∼180 K and are close to that of the H₂O "W" component (T_ex∼200 K; Table 2), which is at −46 km s⁻¹, and we thus conclude that MIR1-W1/MIR2-W2 corresponds to the H₂O "W" component. We emphasize that, because the CO components MIR1-W1 and MIR2-W2 were found to be in front of both binary stars, for the "W" component in H₂O, we do not consider the dilution of the relative intensity (see discussion in Section 2) when performing the rotation diagram analysis in Section 3.

We note that for the "H1" H₂O component, which is centered at −39.5 km s⁻¹ and is close to v_sys, the existence of ν₂ = 2–1 transitions and its derived high temperature exclude a connection between it and the cold CO component ("MIR1-C/MIR2-C"). It is possible that in the low-energy CO transitions, the CO counterpart of this hot H₂O component at −38 km s⁻¹ may be hidden underneath a lower-temperature component at that velocity. We will discuss this possibility in Section 5.2.1.

At high-energy levels, MIR1 and MIR2 contribute differently toward the absorbing components. The component MIR2-H1 in CO is centered at −37.5 km s⁻¹ and is characterized by 600–700 K. We link this to the "H1" water component. It is more difficult to determine the origin of the −54.5 km s⁻¹ "H2" H₂O component. On the one hand, the velocity position, as well as the temperature, of the "H2" H₂O component, are compatible with the CO MIR2-H2 component (Table 4). On the other hand, the "H2" H₂O component is also likely correlated with the hot CO MIR1-W${}^{{\prime} }$ component, which is a complex amalgam of several components. MIR1-W${}^{{\prime} }$ has a complicated origin as is indicated by the varying average line profiles in ¹³CO, ¹²CO, and ¹²CO ν = 2–1. One of the two peaks in the ¹²CO ν = 2–1 profile (at −54 km s⁻¹) of MIR1-W${}^{{\prime} }$ coincides more or less with the −54.5 km s⁻¹ H₂O "H2" component, although the other one (at −39.5 km s⁻¹) has no counterpart in the H₂O spectrum. Likewise, the ¹²CO and ¹³CO MIR1-W${}^{{\prime} }$ components are centered at −46 and −49 km s⁻¹, respectively. Hence, we consider that the H₂O "H2" component is related to one of the CO MIR1-W${}^{{\prime} }$ components.

In conclusion, the water component "H1" is much stronger, and it better matches MIR2-H1 in CO. The water component "H2" is weaker, and its origin is less clear. Thus, we connect "W" to the warm component at −45 km s⁻¹ in CO, "H1" to MIR2-H1, and "H2" to either MIR2-H2 or MIR1-W' . Once the hot components in H₂O are linked to those in CO, we can derive the H₂O/CO abundance ratio (Table 5). Assuming that the CO abundance is equal to the gas-phase C abundance in the diffuse interstellar medium (ISM; 1.6×10⁻⁴; Cardelli et al. 1996; Sofia et al. 1997), we derive H₂O column densities under the disk model of 3.6 × 10¹⁹ and 8.9 × 10¹⁸ cm⁻² for "H1" and "H2," respectively (Table 5).

Table 5. Comparison of H₂O Characteristics Derived from ISO-SWS and SOFIA-EXES Observations

Properties	ISO-SWS a	SOFIA-EXES
		H1 (Slab)	H1 (Disk)	H2 (Slab)	H2 (Disk)
Tex(H2O) (K)	${400}_{-150}^{+200}$	${471}_{-15}^{+14}$	${491}_{-14}^{+13}$	${600}_{-27}^{+28}$	${612}_{-30}^{+27}$
N(H2O) ( cm−2)	${3}_{-1}^{+1}\times {10}^{17}$	${2.5}_{-0.2}^{+0.3}\times {10}^{19}$	3.6 ± 1.2 × 1019	${5.3}_{-0.6}^{+0.3}\times {10}^{18}$	8.9 ± 0.3 × 1018
X[H2O]/X[CO]	0.05	${1.1}_{-0.4}^{+0.4}$	1.5 ± 0.5	${3.8}_{-2.0}^{+1.8}$ or ${0.9}_{-0.4}^{+0.2}$	0.9 ± 0.3

Note.

^a Boonman & van Dishoeck (2003).

Download table as:  ASCIITypeset image

Under the framework of the known structures in W3 IRS 5, the M-band spectroscopic study on CO pictures the kinematic and physical properties of the gaseous environment. Specifically, it includes a shared foreground envelope at −38 km s⁻¹, several high-velocity clumps (referred to as "bullets") from −60 to −100 km s⁻¹, and a few warm and/or hot components in the immediate environment of the binary from −38 to −60 km s⁻¹ (Li et al. 2022). While we have connected "W," "H1," and "H2" in water to MIR1-W1/MIR2-W1, MIR2-H1, and MIR2-H2 (or MIR1-W') in CO (Section 4.2), gaseous water is not detected in the shared cold envelope or the bullets. We refer to past water studies for a more complete view toward W3 IRS 5 (see Table 5) in this section and discuss the implications of these "detections" and "nondetections" from the perspective of both the kinematics and the chemical abundances (Table 6).

5.1. Hot Gaseous Water in W3 IRS 5

5.1.1. Physical Origins of the Hot Components

As we have built the connections of H₂O "H1" and "H2" components to the hot CO components via the kinematic information, it is natural to consider that, as for the corresponding CO components, the two H₂O components have a disk rather than a foreground cloud origin. Although the observed W_λ /λ can be successfully fitted to the theoretical curve-of-growth of either model (see Sections 3 and 4.1), and the two models derive comparable column densities (see Section 4.2) and temperatures, we nevertheless consider the disk model as the more preferred one for the following reasons.

If the hot components have a foreground origin, the heating mechanisms are either due to the radiative heating or due to shocks. On the one hand, if the water components are radiatively heated, the temperature range from 450–600 K corresponds to a distance of 280–140 au in the W3 IRS 5 system (Li et al. 2022). As "H1" is at the systemic velocity, this would imply that this component is static in radial velocity in the highly dynamic environment close to a high-mass protostar. The "H2" component has the opposite issue, as it is moving at a radial velocity of −15 km s⁻¹ relative to the system, and would move outward by some 100 au in 30 yr. As the physical conditions derived from CO observations in 2020 are very similar to those derived in 1991 (Mitchell et al. 1991), this distance range does not work for H2, either, because of its relative velocity of 15 km s⁻¹, and it shall move outward along the line of sight for a large distance (∼100 au). On the other hand, interpreting the warm temperature of these two components as the results of shocks has issues too. For J-shocks, the derived temperatures are quite high. J-shocks initially heat the gas to very high temperatures larger than 10⁴ K. The gas then cools rapidly and, as molecules form, the H₂ formation energy keeps the temperature at ∼400 K for a column density of ∼10²² cm⁻² (Hollenbach et al. 2013). Using the relationship between postshock temperature in this plateau region derived from J shock models (Equation (24) in Hollenbach et al. 2013), a temperature of 500 K (see, Table 5) requires a preshock density of, ${n}_{0}\simeq {10}^{8}{\left(T/500\,{\rm{K}}\right)}^{8.3}\left({v}_{s}/50\,\mathrm{km}/{\rm{s}}\right){\left({\rm{\Delta }}{v}_{D}/1\,\mathrm{km}/{\rm{s}}\right)}^{1.8}$ cm⁻³, where v_s is the shock velocity and Δv_D is the Doppler width in the molecular gas. C-type shocks can yield temperatures of 450–600 K for shock velocity of ∼10 km s⁻¹ (Kaufman & Neufeld 1996). However, a C-shock produces a warm H₂ column density of ∼10²¹ cm⁻², which is far too small to be consistent with the column density we derived for the hot components.

If "H1" and "H2" are in the disks, the H₂O (and CO) excitation temperatures are similar to (but slightly lower than) the dust continuum temperature. Hence, heating is not an issue. However, the problem now resides in the exact positions of the two components on the disk based on their projected velocity information along the line of sight. As pointed out in Li et al. (2022), it is difficult to pinpoint the locations of the blobs on the disk because the inclination angles of both MIR1 and MIR2 and the systematic velocity of MIR1 are unknown. "H1" was connected with "MIR2-H1" because of the similar velocity and temperature. We note that although the central velocities of "H1" and "MIR2-H1" differ by ∼2 km s⁻¹, such a difference is insignificant compared with the uncertainty level of v_LSR of "H1," which is 2.1 km s⁻¹ (see Section 3). For "H2," it is not even clear whether it is associated with MIR1 or MIR2, as the binary protostars are not spatially resolved by SOFIA, and both binary stars show hot components at the same velocity close to −55 km s⁻¹. As shown in Figure 6, in MIR1, the hot component shows up in ¹²CO ν = 2–1 vibrationally excited transitions, which indicate a high-density region, ∼10¹⁰ cm⁻³. In MIR2, the component is possibly a blob on an inclined disk at a distance smaller than 80 au to the central protostar, but such a scenario poses questions on the inclination of the disk MIR2 again. Therefore, we emphasize that disks in MIR1 and MIR2 need to be spatially and spectrally resolved to fully understand the structures in this region. Observations of vibrationally excited lines via submillimeter interferometers may be very instrumental in settling these issues.

Analyses of AFGL 2136 and AFGL 2591 observations have faced the same issues in determining whether the hot 600 K absorbing components reside in the foreground slab or the disk (Barr et al. 2022a). A foreground slab model places very strong constraints on the geometry. In AFGL 2136, if the slab has a water maser origin, the spatial coverage is even larger than the MIR disk, inconsistent with the required covering factor in order to explain the saturation of absorption lines at nonzero flux. In AFGL 2591, wavelength-dependent covering factors are needed to interpret the difference of the spatial coverage derived from the 7 and 13 μm spectroscopy, while the component does not cover the source at all at 3 μm (Barr et al. 2022b). Continuum emission size and chemistry need to be radius dependent, possibly due to the temperature gradient, to interpret the different covering factors. Moreover, similar to W3 IRS 5, too high a density, and too high an abundance argue against a shock origin. In contrast, Keplerian disks as well as clumpy substructures were spatially resolved by Atacama Large Millimeter/submillimeter Array (AFGL 2136; Maud et al. 2019) and NOEMA (AFGL 2591; Suri et al. 2021), supporting the scenario that the absorption arises in blobs in the disk.

In summary, from a broad view, it is a prevalent scenario that hot absorption gas is detected against the MIR continuum backgrounds in massive protostars (e.g., van Dishoeck & Helmich 1996; Cernicharo et al. 1997; Lahuis & van Dishoeck 2000; Boonman & van Dishoeck 2003). Locating a blob on a disk that has a vertical outward-decreasing temperature gradient requires fewer constraints on the geometry than a foreground slab model does. However, disk models face challenges in realizing such an internal heating mechanism. As discussed in Barr et al. (2022a), the flashlight effect may ensure the disk is not externally heated (Nakano 1989; Yorke & Bodenheimer 1999; Kuiper et al. 2010). However, if one proposes the dissipation of gravitational energy, the accretion rate would be orders of magnitude higher than the expected accretion rate (McKee & Tan 2003; Hosokawa et al. 2010; Kuiper et al. 2011; Caratti o Garatti et al 2017). Hence, dissipation of turbulent and/or magnetic energy inherited from the prestellar core would be required, implying a very early and active stage in the formation of these massive protostars.

5.1.2. Vibrationally Excited H₂O

We presented in Figure 3 the rotational temperatures of "H1" and "H2" in the first excited vibrational state (ν₂ = 1), which are 703 ± 60 K and 946 ± 170 K, respectively. Applying curve-of-growth analyses to "H1," the corrected rotational excitation temperatures are ${654}_{-191}^{+135}$ K in the slab model and ${691}_{-212}^{+122}$ K in the disk model (Table 2) and are not far away from the result derived from the rotation diagram. The increment of the total column density of "H1" after correction is ∼3–4 times, much less than that of ∼20 times on the ν₂ = 0 state, indicating a less severe optical depth effect.

One can derive the vibrational excitation temperature T_vib via the Boltzmann equation by comparing the column density in the ν₂ = 0 level, N₀, with N₁ in the ν₂ = 1 level:

$$\begin{eqnarray}&&{N}_{1}/{N}_{0}=\exp (-2294.7\,{\rm{K}}/{T}_{\mathrm{vib}}).\end{eqnarray} \tag{ 16 }$$

We list in Table 7 multiple vibrational excitation temperatures for "H1" and "H2," using N₀ and N₁ before and after corrections on the optical depth effects. Results in Table 7 indicate that such a correction also decreases the derived vibrational excitation temperatures by a few hundred kelvin. Comparing the corrected T_vib with the corrected rotational excitation temperature for the population in the ν₂ = 0 state (Table 2) or in the ν₂ = 1 state, we conclude that those temperatures are in relatively good agreement within the error and that vibrational equilibrium is reached.

Table 7. Vibrational Excitation Temperatures

	H1	H2
Tvib(thin) (K) a	${843}_{-233}^{+553}$	${1127}_{-498}^{+2083}$
Tvib(slab) (K) b	${511}_{-42}^{+30}$	${549}_{-130}^{+86}$
Tvib(disk) (K) c	${581}_{-139}^{+52}$	${488}_{-100}^{+66}$

Notes.

^a Values of N₁ and N₀ in Equation (16) are derived from rotation diagram analysis (see Table 2; same for the columns below). ^b N₁ and N₀ are corrected under the slab model. Since no curve-of-growth analysis is applied for the ν₂ = 1 state of "H2," we adopt N₁ from the rotation diagram analysis, which is the same for the "H2" result under the disk model. ^c N₁ and N₀ are corrected under the disk model.

Download table as:  ASCIITypeset image

The existence of vibrationally excited H₂O implies that the physical conditions of the hot absorbing gas are extreme, as the ν₂ = 1 state lies 2295 K above the ground vibrational state. The two main excitation mechanisms are collisional excitation due to warm, dense gas and radiative excitation by infrared radiation due to warm dust. If the ν₂ = 1 state is collisionally populated, a (postshock) density exceeding 10¹⁰ cm⁻³ is required for thermalization. This order takes account of a critical density ¹⁵ of 10¹¹ cm⁻³ and a radiative trapping effect with β of ∼0.1 for an optically thick line (τ ∼ 10). As a comparison, Barr et al. (2020) estimated a blob density of 10⁹ cm⁻³ in the disk systems of AFGL 2136 and AFGL 2591.

Other than the collisional excitation, one can estimate the relative importance of the excitation due to the strong radiation field. As is described in Tielens (2005), with a dilution factor W (W < 1) on the radiation,

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{ex}}}{{T}_{R}}=1+\displaystyle \frac{{{kT}}_{\mathrm{ex}}}{h\nu }\mathrm{ln}W,\end{eqnarray} \tag{ 17 }$$

in which T_R is the radiation temperature and can be characterized by a dust temperature T_d inside a blackbody. Taking that the hot components are on the surface of the disk and are receiving radiation with W = 0.5 (half of the disk), that T_ex is from 400–500 K (Equation (16)), and that h ν are larger than 2295 K, we derive T_ex/T_R > 0.85. This result indicates that the radiation field does drive the gas to the radiative temperature. We note that at the high implied densities (∼10¹⁰ cm⁻³), collisions between gas and dust will lead to gas kinetic temperatures that are coupled to but slightly lower than the dust temperature (Takahashi et al. 1983).

Li et al. (2022) detected vibrationally excited CO transitions and derived a rotational excitation temperature, 791 K, of CO ν = 1 state for the MIR2-H1 component. The similarity excitation temperatures of the vibrational band of water (703 K) and CO support that the excitation temperature is more representative of the color temperature of the radiation field than the kinetic temperature because water has a much larger dipole moment than CO and therefore much more rapid spontaneous radiation. Finally, we emphasize that the scenario that the kinetic temperature is less than the dust temperature is inherent to our results, as otherwise, one would observe emission rather than absorption lines against the MIR continuum of the observed sources (see Appendix A in Barr et al. 2022a).

5.1.3. Comparison with Past Observations

The temperature range of "H1" and "H2" at −39.5 and −54.5 km s⁻¹ from 400–500 K is comparable with that derived from the ISO-SWS observations in Boonman & van Dishoeck (2003). ¹⁶ In the ISO-SWS study, the individual velocity components were not spectrally resolved, separate transitions were blended (R = 1400, 214 km s⁻¹), and a Doppler width b of 5 km s⁻¹ was assumed in modeling the absorption features. As presented in Table 5, the column density of hot gaseous components derived from our SOFIA/EXES study is about 2 orders of magnitude larger than the ISO-SWS results. This is a significant increment and much more that the 2.4 and 4.3 times increments derived for AFGL 2136 and AFGL 2591 (Barr et al.2022a). We interpret this from two aspects: first, the absorption intensities of W3 IRS 5 measured by ISO-SWS are much lower than those for AFGL 2136 and AFGL 2591 (Boonman & van Dishoeck 2003). This indicates a more significant opacity effect than in AFGL 2136 and AFGL 2591. As a result, the column densities in W3 IRS 5 corrected by the curve-of-growth analysis are much higher than those corrected values in AFGL 2136 and AFGL 2591 (Barr et al. 2022a). Second, saturated lines do not go to 0 but rather reach a nonzero intensity because of either the temperature gradient in a disk atmosphere or a covering factor less than 1 for a foreground cloud.

5.2. Other Foreground Gaseous Components

5.3. Chemical Abundances along the Line of Sight

The availability of data on column densities of different species such as gaseous CO and ices along the line of sight toward W3 IRS 5 makes it an appropriate example to address the oxygen and carbon budget. We discuss below the reservoirs of the two elements in different environments of the W3 IRS 5 system and other massive protostars including the hot disks as well as the cold foreground clouds.

As described in Section 5.1, SOFIA observations derived much higher column densities of hot gaseous water than that of hot gaseous CO in W3 IRS 5 as well as AFGL 2136 and AFGL 2591 (Barr et al. 2022a). While iSHELL measurements (Barr et al. 2020; Li et al. 2022) provide a better constraint on the amount of gaseous CO from the same region, we derive high relative abundances of H₂O to CO; i.e., ∼1–1.5 for W3 IRS 5, 1.6 for AFGL 2136, and 7.4 for AFGL 2591. Such a high relative H₂O to CO abundance is expected for warm, dense gas where gas-phase chemistry rapidly converts the available O not in CO into H₂O (Kaufman & Neufeld 1996). As a comparison, these values are much higher than the H₂O/CO = 10⁻⁴ derived from submillimeter observations by Herschel-HIFI toward the hot core region of AFGL 2591 (Kaźmierczak-Barthel et al. 2014), or the value of 4.4% from the warm 200 K component identified in this water study (see Table 2). On the other hand, a high relative H₂O to CO abundance from 1–2 was observed toward T Tauri and Herbig disks (Carr & Najita 2008; Salyk et al. 2011).

In the cold dense ISM, there is a well-documented problem of the missing oxygen budget (Whittet 2010). While Hollenbach et al. (2009) predicted that oxygen not in silicates or oxides should be eventually converted into gaseous CO and ices, a substantial shortfall of oxygen is observed. In the study of the Taurus complex dark clouds (Whittet 2010), the combined contributions of gaseous CO, ice, and silicate/oxide account for < 300 ppm of the elemental oxygen compared to the solar value of 490 ppm (Asplund et al. 2009). We observed a similar missing oxygen reservoir in W3 IRS 5: assuming 284 ppm of the O in diffuse clouds (Cartledge et al. 2004), we only see a value of 58.1 ppm (20.4%) in cold, dense clouds (see Table 8). If one uses local B stars as the interstellar standard, the total budget of the elemental abundance of oxygen is even higher (575 ppm rather than 490 ppm; Nieva & Przybilla 2012). Therefore, a budget close to the oxygen abundance in silicate (∼200 ppm; Tielens & Allamandola 1987) is missing or is locked up in an unidentified form, which is referred to as "the unidentified depleted oxygen (UDO)" in Whittet (2010).

Table 8. Repository of Elemental Carbon and Oxygen in the Cold Regions of Massive Protostars

	W3 IRS 5			AFGL 2136			AFGL 2591
	N	XC	XO	N	XC	XO	N	XC	XO
	( cm−2)	(ppm)	(ppm)	( cm−2)	(ppm)	(ppm)	( cm−2)	(ppm)	(ppm)
Hydrogen	2.0(23)			7.4(22)			8.0(22)
CO (gas) a	4.7(18)	23.5	23.5	5.1(18)	68.9	68.9	2.5(18)	31.3	31.3
H2O (gas) b	<4.7(15)	⋯	<0.02	⋯	⋯	⋯	⋯	⋯	⋯
CO (ice) c	2.1(17)	1.1	1.1	2.7(17)	3.6	3.6	⋯	⋯	⋯
CO2 (ice) c	7.1(17)	3.6	7.1	7.8(17)	10.5	21.1	1.6(17)	2.0	4.0
H2O (ice) c	5.1(18)	⋯	25.5	5.1(18)	⋯	68.9	1.2(18)	⋯	15.0
CH3OH (ice) c	1.7(17)	0.9	0.9	2.6(17)	3.5	3.5	1.7(17)	2.1	2.1
Sum (ice)	⋯	5.6	34.6	⋯	17.6	97.1	⋯	4.1	21.1
Sum (total)	⋯	29.1	58.1	⋯	86.5	166.0	⋯	35.4	52.4
Diffuse Clouds (gas)		160 d	284 e		160	284		160	284
Silicate Abundance f			200			200			200
Solar Abundance g		269	490		269	490		269	490
B-stars Abundance h		214	575		214	575		214	575

Notes. (1) For column densities, powers of 10 are given in parentheses. (2) X_C and X_O are the relative abundances derived from N_C/N_H and N_O/N_H.

^a Results of gaseous CO are adopted from Li et al. (2022) for W3 IRS 5, and are from Barr et al. (2020) for AFGL 2136 and AFGL 2591. The temperatures of the cold regions are ∼50 K, 27 K, and 49 K for W3 IRS 5, AFGL 2136, AFGL 2591, separately (Barr et al. 2020; Li et al. 2022). ^b Results from this work. ^c All measurements on ice are adopted from Gibb et al. (2004). ^d Cardelli et al. (1996), Sofia et al. (1997). ^e Cartledge et al. (2004). ^f Tielens & Allamandola (1987). ^g log_C = 8.43, log_O = 8.69 (Asplund et al. 2009). ^h Nieva & Przybilla (2012).

Download table as:  ASCIITypeset image

The intrinsic properties of the reservoir of the missing oxygen remain mysterious. Refractory dust compounds like carbonates are implausible, as they will also survive and appear in the diffuse medium. Neither gas-phase nor solid-phase O₂ are possible as well. While SWAS observations (Goldsmith et al. 2000) toward massive protostars provide an upper limit of 0.1 ppm on the gaseous O₂, solid O₂ is too volatile in the line of sight of W3 IRS 5. Oxygen-bearing organics in solid ice would be a potential carrier, although no significant detection of the related bonds has yet been detected at infrared wavelengths (Gibb et al. 2004). We suggest that the MIRI spectrograph on board the JWST is well suited to study such an organic inventory. One other possibility of the UDO is a population of very large water ice grains (>1 μm) in the cold gas (Jenkins 2009). These large grains produce gray extinction, and their putative presence is hard to refute.

Similar to the "oxygen crisis," a depletion problem in elemental carbon exists in the envelope of W3 IRS 5 (see Table 8). The total amount of carbon (32 ppm) comprises only 19.7% of the value expected in diffuse clouds (160 ppm; Cardelli et al. 1996; Sofia et al. 1997). As discussed in Li et al. (2022), one speculation is that the carbon-containing ice compounds were converted into an organic residue by prolonged UV photolysis (Bernstein et al. 1995, 1997; Vinogradoff et al. 2013).

Both the "oxygen crisis" and the "carbon crisis" were observed in AFGL 2136 and AFGL 2591 as well. In contrast to previous studies that relied on comparing pencil beam IR absorption line studies with submillimeter emission observations, the IR pencil beam samples the same material in absorption for these massive protostars. As shown in Table 8, the depletion problems are more severe for AFGL 2591, but less severe for AFGL 2136. We suggest that further studies of the different oxygen reservoirs could help pinpoint the processes involved in the missing oxygen or carbon reservoirs by studying a large enough sample with diverse characteristics.

We conducted high spectral resolution (R = 50,000; 6 km s⁻¹) spectroscopy from 5–8 μm with EXES on board SOFIA toward the hot core region associated with the massive binary protostar W3 IRS 5. By comparing with the LTE models constructed with the existing laboratory line information, we identified about 180 ν₂ = 1 − 0 and 90 ν₂ = 2 − 1 absorption lines. Preliminary Gaussian fittings and rotation diagram analyses reveal two hot components with T > 600 K and one warm component of 190 K. However, the large scatter in the rotation diagrams of the two hot components reveals (1) opacity effects, (2) that the absorption lines are not optically thin, and (3) the total column densities derived from the rotation diagrams are underestimated.

We adopted two curve-of-growth analyses to account for the opacity effects of the hot components. One model considers absorption in a foreground slab partially covering the background emission. The other model assumes absorption in the photosphere of a circumstellar disk with an outward-decreasing temperature in the vertical direction. In both models, about half of the data points converted from the ν₂ = 1 − 0 transitions are located on the logarithmic part of the curve-of-growth, confirming that the corresponding absorption lines are optically thick. The two curve-of-growth analyses correct the column densities by at least an order of magnitude and lower the derived excitation temperatures accordingly. We note that for the disk model, the results of the curve-of-growth analysis depend on the adopted velocity width σ_v and the parameter that characterizes absorption and scattering of the absorption line. We provide a reference table for different σ_v and , as the two parameters are poorly constrained.

Although our SOFIA-EXES observations do not spatially resolve the binary protostars in W3 IRS 5, using the kinematic and temperature characteristics, we link each H₂O component to a spatially separate CO component identified in IRTF/iSHELL observations (R = 88,100). Specifically, the warm H₂O component "W" is linked to the shared warm CO component MIR1-W1/MIR1-W2, the hot H₂O component "H1" is linked to the MIR2-H2 in CO, and the hot H₂O component "H2" is considered to be related to one of the CO MIR1-W${}^{{\prime} }$ components.

Once the connections of H₂O components and CO components were established, we discussed the physical origins of H₂O components in light of a better understanding of CO components. From our analysis, we conclude that the disk model is the preferred one over the slab model, due to consideration of the geometry constraints, although disk models face challenges in realizing such an internal heating mechanism.

We derive the H₂O/CO abundance ratio based on the results of the disk model and discuss the chemical abundances along the line of sight based on the H₂O-to-CO connection. For the hot gas, we derive a high H₂O/CO abundance ratio of 0.9. Such a high relative H₂O to CO abundance is expected for warm, dense gas where gas-phase chemistry rapidly converts the available O not in CO into H₂O. For the cold gas, we observe a substantial shortfall of oxygen in agreement with earlier studies of cold dense clouds. We suggest that organics in solid ice are the potential carrier.

We thank the anonymous referee for providing thorough comments and helpful suggestions that greatly improved the quality of this paper. We acknowledge the support for the EXES Survey of the Molecular Inventory of Hot Cores (SOFIA No. 08-0136) at the University of Maryland from NASA (NNA17BF53C) Cycle Eight GO Proposal for the Stratospheric Observatory for Infrared Astronomy (SOFIA) project issued by USRA.

Software: Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020).

We present in Figure 7 the important telluric lines from 5.36–7.92 μm produced by the PSG models under representative observational parameters (see Table 9). Eight molecular species, H₂O, O₃, CH₄, NO₂, N₂O, CO₂, HNO₃, and O₂, are among the most important ones. For the 15 observational settings, the input surface pressure and scaling factors are tuned by hand to visually match the telluric features in observed spectra, and therefore, may not represent the actual values.

Zoom In Zoom Out Reset image size

The standard star, Sirius, was observed only in one setting from 7.28–7.46 μm. While the Sirius spectrum has advantages in reflecting the actual baseline, we compare the qualities of the results derived by using the Sirius spectrum as well as using the PSG models. The data reduction with Sirius is straightforward, as the median filtering processes or the modeling of telluric lines are not needed.

As a result, we present in Figure 8 the equivalent widths of each individually identified line derived from the two data reduction methods. We conclude that the results are in good agreement with each other by ∼10%.

Zoom In Zoom Out Reset image size

As presented in Sections 3 and 4, data analyses in this paper are based on the decomposition of the absorption profiles into the three components at −54.5, − 45, and −39.5 km s⁻¹ in low-energy levels, while the −39.5 km s⁻¹ component is hot. However, we realize that one cold component of ∼50 K may exist at −38 km s⁻¹ because of the detection of the cold CO component (Li et al. 2022).

We argue that a cold component is possibly hidden in Figure 9, in which a list of accumulative average spectra is presented. We note that the absorption feature at −35 km s⁻¹ is possibly related to the cold component because it disappears as the energy levels increase. We also present two different Gaussian fitting methods for the average spectrum below 200 K in Figure 10, while in one, we fix the right wing with a component at −38 km s⁻¹ and in the other at −40 km s⁻¹. We conclude that the previous one provides a better fit, suggesting that at this energy level, the cold −38 km s⁻¹ component possibly dominates the line profile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even if the cold −38 km s⁻¹ component does exist, we conclude that the water-to-CO abundance is still low at this temperature. Otherwise, very saturated water absorption lines will dominate the line profiles. We estimate an upper limit of such a water-to-CO ratio of 0.4% by Figure 11.

Zoom In Zoom Out Reset image size

Figure 12 presents the two grid-search results and the best-fitted curve-of-growth for the "H2" component. Figure 13 presents the results for the ν₂ = 2–1 transitions on the "H1" component.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We present from Table 10–14 the properties of decomposed water lines, H₂O ν₂ = 1–0 and ν₂ = 2–1, of W3 IRS 5 from 5.36–7.92 μm in this study. For "H2" in the ν₂ = 1–0 transitions, the velocity centers of lines with E_l < 800 K are fixed at −54 km s⁻¹. For "W," in the ν₂ = 1–0 transitions, the velocity centers are fixed at −45 km s⁻¹. τ_p,thin, τ_p,slab, and β₀ are calculated from the best-fitted temperatures and column densities from the rotation diagram, the slab model, and the disk model.