The H₂O Spectrum of the Massive Protostar AFGL 2136 IRS 1 from 2 to 13 μm at High Resolution: Probing the Circumstellar Disk

CRIRES observations targeting the ν₁ and ν₃ rovibrational bands of H₂O were made at two reference wavelengths: 2480.0 nm and 2502.8 nm. Due to the lack of a bright natural guide star near AFGL 2136, the adaptive optics system was not utilized. The 02 slit was employed, providing a spectral resolving power of about 10⁵, corresponding to ∼3 km s⁻¹ resolution. The slit was oriented at a position angle of 45° (along a northeast-to-southwest axis) to minimize contamination from the near-infrared scattered light observed in the region (e.g., Murakawa et al. 2008). Spectra were obtained in an ABBA pattern along the slit with 10'' between the two nod positions and ±3'' jitter width. All observations of AFGL 2136 IRS 1 were immediately preceded by observations of the bright A2IV star HR 6378 using the same setup for use as a telluric standard.

EXES observations targeting the ν₂ rovibrational band of H₂O were made in cross-dispersed high-resolution mode targeting central wavenumbers of 1485.24 cm⁻¹, 1639.29 cm⁻¹, and 1747.25 cm⁻¹ (hereafter referred to as the 6.7 μm, 6.1 μm, and 5.7 μm spectra, respectively). The entrance slit had a width of 184, providing a resolving power (resolution) of ∼85,000 (3.5 km s⁻¹), and a length of about 10'' (varies slightly between settings). To facilitate the removal of telluric emission lines, exposures alternated between on-target and a blank sky position 15'' away. AFGL 2136 was one of six science targets observed as part of SOFIA program 04-0120 over multiple flights. The bright A1V star α CMa (Sirius) was observed only once at the same three settings given above for use as a telluric standard for our entire program, so it was not necessarily observed on the same night as AFGL 2136 IRS 1. As will be shown here—and in the subsequent paper discussing the other sources—this strategy proved sufficient for removing telluric features from the spectra.

TEXES observations targeting pure rotational transitions of H₂O (and rovibrational transitions of other molecular species) were made in high-med mode (cross-dispersed echelon + 32 l/mm echelle) at reference wavenumbers of 931.5 cm⁻¹, 856.0 cm⁻¹, 806.5 cm⁻¹, and 768.5 cm⁻¹. A resolving power (resolution) of ∼85,000 (3.5 km s⁻¹) was achieved using the 054 wide slit. The target was nodded 16 within the ∼4'' long slit between exposures to facilitate the removal of sky emission lines. Immediately before or after observations of AFGL 2136 IRS 1, the asteroid 16 Psyche was observed using the same configurations as above for use as a telluric standard.

The following reduction procedure was performed independently for each of the four CRIRES detectors on each of the six nights that observations occurred. Raw images were processed using the CRIRES pipeline version 2.3.3. Standard calibration techniques, including subtraction of dark frames, division by flat fields, interpolation over bad pixels, and correction for detector nonlinearity effects, were applied. Consecutive A and B nod position images were subtracted from each other to remove sky emission features, and all images from each nod position were combined to create average A and B images. Spectra were extracted from these images using the apall routine in iraf⁸ and then imported to IGOR Pro.⁹ Wavelength calibration was performed using atmospheric absorption lines, and is accurate to ±1 km s⁻¹. Spectra from the A and B nod positions were then averaged onto a common wavelength scale, and were further processed using the shared reduction steps described in Section 3.4.

Data were processed using the Redux pipeline (Clarke et al. 2015) with the fspextool software package—a modification of the Spextool package (Cushing et al. 2004)—which performs source profile construction, extraction and background aperture definition, optimal extraction, and wavelength calibration for EXES data. We used this software to produce wavelength-calibrated spectra for each individual order of the echellogram. These individual spectra were then stitched together using an average of both orders in the overlap regions to produce a continuous spectrum for each of the three separate observations. Further processing is described in Section 3.4.

Data were processed using the TEXES pipeline (Lacy et al. 2002). The pipeline performs standard calibration techniques, including flat-field calibration, removal of signal spikes, interpolation over dead and/or noisy pixels, subtraction of nod pair images, summation of differenced spectral images, spectral extraction, and wavelength calibration. Output products are wavelength-calibrated spectra for each echelle order. These individual echelle order spectra were then combined onto a common wavelength scale to create a single spectrum for each setting. The 856.0 cm⁻¹, 806.5 cm⁻¹, and 768.5 cm⁻¹ settings have gaps in wavelength coverage between adjacent echelle orders, so the intensities from the individual orders are unchanged by this combination. The 931.5 cm⁻¹ setting produces echelle orders that overlap in wavelength, so an average of both orders in the overlap regions was again used to produce a continuous spectrum. Further processing is described in Section 3.4.

Once data were processed to the point where each entry in Table 1 had a single, wavelength-calibrated spectrum, subsequent data reduction methods for all spectra were shared, regardless of the instrument used. To remove baseline fluctuations and atmospheric features, the science target spectra were divided by the corresponding telluric standard spectra using custom macros developed in IGOR Pro that allow for stretching and shifting of the telluric standard spectrum in the wavelength axis, and scaling of the telluric standard intensity according to Beer's law (McCall 2001). The resulting ratioed spectra were then divided by a 30-pixel boxcar average of the continuum level (interpolated across absorption lines) to remove residual fluctuations and produce normalized spectra. Regions where strong atmospheric absorption features result in no useful information were removed from the spectra to improve visualization. Because Earth's orbital motion causes astrophysical lines to shift with respect to atmospheric lines (i.e., observed wavelength) throughout the year, wavelength scales for all spectra were converted to the local standard of rest (LSR) frame.

One final step was required for processing the CRIRES observations since data were acquired at the same settings on multiple nights. Where multiple spectra cover the same wavelength range, they were combined using a weighted average (weighted by 1/σ², where σ is the standard deviation of the line-free continuum in each spectrum). In this way, all of the individual CRIRES spectra were combined to form a single normalized spectrum. This spectrum and the normalized spectra resulting from the EXES and TEXES observations are presented in Figures 1–3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

More than 100 absorption features are present in these spectra, including many that are blends due to absorption from multiple transitions. Using the temperature and total H₂O column density inferred from our previous study of this source (Indriolo et al. 2013a) and assuming local thermodynamic equilibrium (LTE), we predicted the strength of all H₂O absorption lines covered by our data. By comparing these predictions to the observed spectra, we identified absorption features that are due to individual transitions of H₂O (i.e., are not blends of multiple H₂O transitions or blends of an H₂O transition and that of another species). Starting with some of the stronger ν₃ band features, we fit absorption lines (in optical depth space) with a sum of two Gaussian components to determine average fit parameters. The line-center velocities (v_cent) and line widths (σ_v) determined from these fits were then used as initial guesses for fitting the other unblended absorption features. During the fitting procedure, these free parameters were restricted to values similar to the initial guesses. The line-center velocities of both components were constrained to remain within ±4 km s⁻¹ of the initial guesses (20–28 km s⁻¹ and 28–36 km s⁻¹ for components 1 and 2, respectively), and line widths for both components were constrained to the range 1 km s⁻¹ ≤ σ_v ≤ 6 km s⁻¹. In this manner, all of the unblended H₂O absorption features were fit with line profiles that are the sum of two Gaussian functions in optical depth. Some example fits are shown in Figure 4. Note that it is unlikely that the two components utilized in the fit correspond to two distinct physical components, as will be discussed in Section 5.4.

Zoom In Zoom Out Reset image size

Using the best-fit parameters, we determined an integrated optical depth ($\int \tau {dv}$) for each component. Uncertainties were computed using the covariance matrix returned by scipy.optimize.curve_fit. The column density in the lower rotational state associated with the transition (N_l) was calculated via the standard equation for optically thin absorption,

$$\begin{eqnarray}&&{N}_{l}=\displaystyle \frac{8\pi }{{A}_{{ul}}{\lambda }^{3}}\displaystyle \frac{{g}_{l}}{{g}_{u}}\int \tau {dv},\end{eqnarray} \tag{ 1 }$$

where λ is the transition wavelength, A_ul is the transition spontaneous emission coefficient, and g_l and g_u are the statistical weights of the lower and upper states, respectively, including nuclear spin degeneracy (i.e., g = (2J + 1)(2I + 1), where I = 0 for para-H₂O and I = 1 for ortho-H₂O). Transition data were taken from the HITRAN database (Rothman et al. 2013). Total column densities in each rotational state are taken to be the sum of the column densities determined from the two components, and the resulting $\mathrm{ln}({N}_{l}/{g}_{l})$ are presented in Table 2.

Table 2.  Results from Unblended H₂O Absorption Lines

v'−v''	${J}_{{K}_{a},{K}_{c}}^{\prime} -{J}_{{K}_{a},{K}_{c}}^{\prime\prime}$	Wavelength	El	Aul	gu	gl	$\int \tau {dv}$	ln $({N}_{l}/{g}_{l})$	${N}_{\mathrm{meas}}/{N}_{\mathrm{pred}}$
		(μm)	(K)	(s−1)			(km s−1)
ν3 = 1 − 0	83,5–71,6	2.446108	1013.2	0.979	51	45	0.79 ± 0.15	35.53 ± 0.19	0.96
ν3 = 1 − 0	55,1–43,2	2.447276	550.4	0.082	33	27	0.18 ± 0.16	36.95 ± 0.90	1.58
ν3 = 1 − 0	72,5–60,6	2.450228	642.7	0.321	15	13	0.20 ± 0.15	36.52 ± 0.75	1.23
ν3 = 1 − 0	84,4–72,5	2.451345	1125.7	1.480	51	45	0.87 ± 0.15	35.22 ± 0.17	0.88
ν3 = 1 − 0	145,9–135,8	2.453230	3783.0	20.510	87	81	0.19 ± 0.17	30.50 ± 0.92	1.56
ν3 = 1 − 0	64,3–52,4	2.454677	598.8	0.339	13	11	0.13 ± 0.14	36.16 ± 1.04	0.78
ν3 = 1 − 0	144,10–134,9	2.461582	3645.6	25.020	87	81	0.22 ± 0.16	30.46 ± 0.72	1.13
ν1 = 1 − 0	84,5–71,6	2.467807	1013.2	0.194	51	45	0.27 ± 0.16	36.07 ± 0.58	1.64
ν1 = 1 − 0	126,7–115,6	2.468616	2876.1	6.393	75	69	0.22 ± 0.17	31.96 ± 0.76	1.10
ν3 = 1 − 0	134,9–124,8	2.468911	3173.4	25.930	27	25	0.15 ± 0.15	31.24 ± 0.98	0.97
ν3 = 1 − 0	73,4–61,5	2.472935	781.1	1.202	15	13	0.51 ± 0.15	36.09 ± 0.30	1.05
ν1 = 1 − 0	116,5–105,6	2.473202	2472.8	2.726	69	63	0.18 ± 0.16	32.69 ± 0.87	1.02
ν3 = 1 − 0	146,8–136,7	2.473837	3965.9	16.270	87	81	0.11 ± 0.17	30.16 ± 1.62	1.60
ν1 = 1 − 0	115,6–104,7	2.475249	2275.2	2.115	69	63	0.22 ± 0.20	33.17 ± 0.89	1.12
ν3 = 1 − 0	143,11–133,10	2.477481	3474.2	24.950	87	81	0.32 ± 0.16	30.82 ± 0.51	1.16
ν3 = 1 − 0	124,8–114,7	2.478598	2732.2	26.170	75	69	0.95 ± 0.15	32.01 ± 0.16	0.87
ν1 = 1 − 0	93,6–82,7	2.479768	1274.2	0.170	57	51	0.17 ± 0.17	35.57 ± 1.05	1.68
ν1 = 1 − 0	54,1–41,4	2.480526	323.5	0.049	33	27	0.14 ± 0.17	37.17 ± 1.23	1.24
ν3 = 1 − 0	133,10–123,9	2.483376	3029.9	25.500	27	25	0.20 ± 0.15	31.48 ± 0.79	0.92
ν3 = 1 − 0	62,4–50,5	2.484255	468.1	0.494	39	33	1.05 ± 0.14	36.73 ± 0.14	1.08
ν1 = 1 − 0	106,4–95,5	2.486180	2122.2	2.179	21	19	0.12 ± 0.16	33.66 ± 1.40	1.34
ν2 = 2 − 0	153,12–142,13	2.488795	3349.3	8.645	93	87	0.18 ± 0.17	31.23 ± 0.91	1.37
ν3 = 1 − 0	114,7–104,6	2.490420	2325.7	24.690	23	21	0.70 ± 0.17	32.93 ± 0.24	0.97
ν3 = 1 − 0	134,10–124,9	2.491232	3057.3	24.010	81	75	0.55 ± 0.16	31.46 ± 0.28	0.96
ν3 = 1 − 0	44,0–32,1	2.492091	305.2	0.127	27	21	0.27 ± 0.14	37.10 ± 0.51	1.12
ν3 = 1 − 0	63,3–51,4	2.494654	574.7	0.922	39	33	1.27 ± 0.14	36.28 ± 0.11	0.85
ν3 = 1 − 0	113,8–103,7	2.496575	2213.0	27.470	23	21	0.67 ± 0.15	32.77 ± 0.22	0.66
ν3 = 1 − 0	132,11–122,10	2.497063	2820.3	27.190	27	25	0.22 ± 0.13	31.52 ± 0.61	0.63
ν3 = 1 − 0	122,10–112,9	2.504591	2432.5	27.810	75	69	1.60 ± 0.15	32.44 ± 0.09	0.73
ν3 = 1 − 0	114,8–104,7	2.511261	2275.2	23.560	69	63	1.81 ± 0.15	32.80 ± 0.08	0.77
ν3 = 1 − 0	112,9–102,8	2.512036	2068.9	27.820	23	21	1.17 ± 0.15	33.30 ± 0.13	0.84
ν3 = 1 − 0	115,7–105,6	2.512708	2472.8	19.930	69	63	1.27 ± 0.16	32.62 ± 0.12	0.95
ν3 = 1 − 0	93,6–83,5	2.517545	1510.9	28.400	19	17	2.33 ± 0.14	34.15 ± 0.06	0.65
ν3 = 1 − 0	116,6–106,5	2.518460	2697.7	17.010	69	63	0.46 ± 0.15	31.76 ± 0.33	0.63
ν3 = 1 − 0	102,8–92,7	2.519438	1729.3	28.430	63	57	3.95 ± 0.13	33.48 ± 0.03	0.51
ν2 = 1 − 0	74,4–73,5	5.705113	1175.0	2.981	15	15	3.97 ± 0.50	34.72 ± 0.12	0.59
ν2 = 1 − 0	81,8–70,7	5.709642	843.5	11.840	51	45	15.25 ± 1.18	33.46 ± 0.08	0.09
ν2 = 1 − 0	113,8–112,9	5.713061	2432.5	3.049	69	69	1.71 ± 0.47	32.33 ± 0.27	0.66
ν2 = 1 − 0	64,3–63,4	5.716237	933.7	2.772	39	39	8.11 ± 0.52	34.55 ± 0.06	0.31
ν2 = 1 − 0	52,4–41,3	5.718679	396.4	5.157	11	9	9.94 ± 6.93	35.39 ± 0.70	0.24
ν2 = 2 − 1	101,10–90,9	5.720222	3614.6	25.810	63	57	2.04 ± 0.47	30.46 ± 0.23	1.06
ν2 = 1 − 0	54,2–53,3	5.721717	725.1	2.390	11	11	5.76 ± 0.50	35.62 ± 0.09	0.59
ν2 = 1 − 0	81,7–80,8	5.726424	1070.5	1.397	17	17	2.48 ± 0.46	34.87 ± 0.19	0.56
ν2 = 1 − 0	44,0–43,1	5.728111	552.3	1.617	9	9	4.28 ± 0.48	35.91 ± 0.11	0.56
ν2 = 1 − 0	83,6–82,7	5.731998	1274.2	2.692	51	51	6.38 ± 0.53	34.06 ± 0.08	0.37
ν2 = 1 − 0	32,1–31,2	6.075447	249.4	6.715	21	21	12.62 ± 0.50	34.54 ± 0.04	0.08
ν2 = 2 − 1	51,4–50,5	6.080180	2763.5	6.406	33	33	1.78 ± 0.58	32.18 ± 0.33	1.09
ν2 = 1 − 0	53,2–44,1	6.088700	702.3	0.351	33	27	3.48 ± 0.40	35.75 ± 0.12	0.64
ν2 = 1 − 0	52,4–43,1	6.096408	552.3	0.733	11	9	3.33 ± 0.34	36.06 ± 0.10	0.65
ν2 = 2 − 1	22,1–21,2	6.100969	2412.9	5.916	15	15	0.84 ± 0.32	32.28 ± 0.38	0.60
ν2 = 1 − 0	53,3–44,0	6.103987	702.3	0.343	11	9	1.72 ± 0.38	36.16 ± 0.22	0.96
ν2 = 1 − 0	42,2–33,1	6.106193	410.4	0.529	9	7	2.73 ± 0.31	36.38 ± 0.11	0.67
ν2 = 1 − 0	31,2–22,1	6.106826	194.1	1.123	21	15	9.60 ± 0.72	36.04 ± 0.08	0.31
ν2 = 2 − 1	40,4–31,3	6.113164	2502.7	15.800	9	7	1.44 ± 0.36	32.34 ± 0.25	0.77
ν2 = 1 − 0	11,1–00,0	6.116331	0.0	7.457	3	1	12.55 ± 0.66	36.36 ± 0.05	0.29
ν2 = 1 − 0	103,8–104,7	6.705153	2275.2	9.459	63	63	4.03 ± 0.27	31.66 ± 0.07	0.25
ν2 = 1 − 0	40,4–51,5	6.707693	469.9	8.449	9	11	5.46 ± 0.24	34.03 ± 0.04	0.07
ν2 = 2 − 1	42,2–43,1	6.711310	2886.1	14.320	9	9	0.89 ± 0.25	31.69 ± 0.28	0.85
ν2 = 1 − 0	61,6–62,5	6.712122	795.5	4.825	39	39	6.10 ± 0.33	33.23 ± 0.05	0.06
ν2 = 1 − 0	63,4–72,5	6.714553	1125.7	1.928	39	45	5.91 ± 0.26	34.11 ± 0.04	0.29
ν2 = 1 − 0	82,7–83,6	6.715692	1447.6	7.602	51	51	6.52 ± 0.25	32.57 ± 0.04	0.12
ν2 = 1 − 0	134,10–135,9	6.722649	3721.4	11.270	27	27	0.41 ± 0.23	30.04 ± 0.56	0.87
ν2 = 1 − 0	21,1–32,2	6.723374	296.8	8.223	5	7	5.22 ± 0.58	34.59 ± 0.11	0.09
ν2 = 2 − 1	71,6–72,5	6.726109	3442.4	20.140	45	45	1.59 ± 0.27	30.30 ± 0.17	0.65
ν2 = 2 − 1	32,1–33,0	6.733401	2744.6	9.039	21	21	1.21 ± 0.26	31.59 ± 0.21	0.59
ν2 = 1 − 0	105,6–114,7	6.737377	2732.2	0.804	63	69	0.34 ± 0.23	31.65 ± 0.67	0.60
ν2 = 2 − 1	73,4–74,3	6.753573	3700.7	17.650	45	45	1.07 ± 0.26	30.03 ± 0.24	0.82
	163,14–150,15	11.724550	3393.1	1.682	99	93	0.26 ± 0.03	28.52 ± 0.13	0.10
	164,13–151,14	12.375657	3785.8	4.211	99	93	0.99 ± 0.05	28.78 ± 0.05	0.28
	163,13–152,14	12.407046	3785.8	4.217	33	31	0.35 ± 0.06	28.82 ± 0.17	0.29
	165,12–152,13	13.033327	4132.6	6.761	99	93	0.64 ± 0.06	27.71 ± 0.09	0.19
H218O
ν2 = 1 − 0	21,2–10,1	6.074293	34.2	6.533	15	9	0.84 ± 0.34	32.20 ± 0.40	...
ν2 = 1 − 0	32,1–31,2	6.103487	248.7	6.404	21	21	1.29 ± 0.50	32.30 ± 0.39	...

Note. Transition properties for all unblended lines fit during our analysis are shown here. Integrated optical depths ($\int \tau {dv}$) are over the entire absorption feature (i.e., the sum of the two components). Individual state column densities assuming optically thin absorption are presented as ln(N_l/g_l). These values should be considered with extreme caution, given the discussion presented in this paper. The last column gives the ratio between the measured and predicted state-specific column densities. Predicted column densities assume an isothermal absorbing slab in LTE with the temperature and total H₂O column density determined from CRIRES observations alone (see Table 3). Smaller values correspond to more "underpopulated" states.

Download table as:  ASCIITypeset images: 1 2

Column densities determined for individual rotational states are used to construct a rotation (Boltzmann) diagram for H₂O (Figures 5 and 6). For gas in LTE at a single temperature (T), the points should form a straight line with a slope equal to −T⁻¹. While this is reasonably true in Figure 5 panel (a), the large scatter in $\mathrm{ln}({N}_{l}/{g}_{l})$ is not explained by this simple isothermal absorbing slab model. In particular, the column densities of a significant number of states fall below the best-fit straight line for the majority of levels by as much as a factor of 10. To investigate this issue, we explore whether or not certain lower-state and/or transition characteristics correlate with states that appear to be "underpopulated."

Zoom In Zoom Out Reset image size

In Figure 5, the points on the rotation diagram are distinguished by different discrete properties of the lower rotational state or the observed transition. Panel (b) has points separated by the instrument with which each transition was observed. It is clear that points from different instruments behave in different ways, but it is unlikely that this is an effect of the different instruments themselves (see Section 5.3). The SOFIA/EXES data show significant scatter, so that subset of data can be used to further investigate the observed behavior. Panel (c) shows EXES data separated by the change in the rotational quantum number of the transition (ΔJ = 1, 0, and −1 are the R, Q, and P branches, respectively). There does not appear to be any significant correlation between ΔJ and $\mathrm{ln}({N}_{l}/{g}_{l})$. Panel (d) shows EXES data separated by the nuclear spin (the ortho and para nuclear spin configurations correspond to ${K}_{a}+{K}_{c}$ being odd and even, respectively) of the lower state, and again, there is no correlation with $\mathrm{ln}({N}_{l}/{g}_{l})$. Although not presented here, the CRIRES data similarly show no correlation between $\mathrm{ln}({N}_{l}/{g}_{l})$ and nuclear spin configuration, and indicate an ortho-to-para ratio of 3:1. Panel (e) separates the EXES data by the three different observed wavelength ranges. While there is some distinction between the different points, "underpopulated" states are not unique to any given wavelength range, and above about 2500 K there is no discrepancy between any of the points. As such, the scatter does not appear to be caused by the three separate EXES observations at different wavelengths.

In Figure 6, the points on the rotation diagram are distinguished by color scales that are used to indicate the spontaneous emission coefficient of the transition from which each column density was determined (panels (a), (c), and (e)) and the transmission level at the point of deepest absorption within each absorption feature from which column densities were determined (panels (b), (d), and (f)). For EXES observations, the "underpopulated" states have column densities determined from stronger absorption features (panel (b)), and for states with E_l < 2500 K, the scatter in $\mathrm{ln}({N}_{l}/{g}_{l})$ is correlated with A_ul, with column densities in "underpopulated" states derived from transitions with larger spontaneous emission coefficients (panel (a)). Hints of these same trends are possibly present in the CRIRES observations shown in panels (c) and (d), although they are within the range of uncertainties for many of the inferred column densities. Panels (e) and (f) show the same plots for all data. It is clear from panel (f) that there is a correlation between the depth of absorption features and how far the column densities derived from those features deviate from their expected values. Only the TEXES data points deviate from this pattern. A more in-depth discussion about what may cause these behaviors follows in Section 5.3.

Zoom In Zoom Out Reset image size

Despite the scatter present in the rotation diagrams, a linear fit to the data points can still be used to infer the gas temperature and total H₂O column density. It is clear from Figure 5, though, that this fit will change based on the subset of data points being used, as is demonstrated by the results presented in Table 3. Given the large scatter in the EXES data and the relatively linear, self-consistent behavior of the CRIRES data, we adopt the results obtained from the CRIRES data alone: $N({{\rm{H}}}_{2}{\rm{O}})=(8.25\pm 0.95)\times {10}^{18}$ cm⁻² and T = 502 ± 12 K. A fit using the 6.7 μm EXES data alone results in a temperature that is 1.6 times higher than, and a total H₂O column density that is 8 times lower than, the adopted values. This implies that results derived from different and/or limited wavelength ranges should be viewed with caution (e.g., Indriolo et al. 2015), and that future observations should be planned with care to minimize such effects. Table 3 also shows that a fit to all of the EXES data does not agree with our adopted values. As the upper envelope of EXES data points is consistent with the CRIRES data though (panel (b) of Figure 5), it should still be possible to infer reasonably accurate values of T and N(H₂O) in sources where only EXES observations have been made by restricting the fit to a subset of the EXES data points. For example, limiting the fit to EXES data where column densities were derived from weak absorption features (line-center transmission >75%; τ₀ < 0.29) produces results in better agreement with the CRIRES results (see Table 3). In future cases where scatter similar to that shown here is seen in rotation diagrams, we recommend restricting the linear fit to only include state-specific column densities determined from weak absorption lines when deriving total column densities and rotation temperatures.

4.3.1. Absorption Lines

A variety of molecular absorption features are observed toward AFGL 2136 IRS 1. Some of these species (NH₃, HCN, and C₂H₂) dominate the absorption features at 10–13 μm, but at 2.5 μm, H₂O is the most prevalent absorber. Despite the fact that a uniform absorbing slab model cannot successfully predict the H₂O absorption from 2 to 13 μm, such a model is still useful in analyzing the CRIRES observations to search for features caused by other species. Using the average Gaussian parameters inferred from the two-component fit described previously (v₁ = 24.8 km s⁻¹, σ₁ = 4.16 km s⁻¹, v₂ = 33.7 km s⁻¹, σ₂ = 2.92 km s⁻¹), adopting the total $N({{\rm{H}}}_{2}{\rm{O}})$ and T inferred from only the CRIRES data, and assuming that the component at 24.8 km s⁻¹ contains 70% of the material, we generate a synthetic H₂O absorption spectrum (see Figures 1, 7, and 8). The top panel of Figure 7 shows the synthetic H₂O spectrum at the wavelengths of the v = 1–0 R(0), R(1), and R(2) transitions of HF, demonstrating that all three HF transitions are detected in absorption. The bottom panel of Figure 7 shows the AFGL 2136 IRS 1 spectrum after division by the synthetic H₂O spectrum, effectively showing the absorption due to only HF. By fitting the HF absorption features after division by the synthetic H₂O spectrum, we find individual state column densities of N(0) = (3.1 ± 0.8) × 10¹⁴ cm⁻², N(1) = (2.6 ± 1.3) × 10¹⁴ cm⁻², and N(2) = (1.6 ± 1.4) × 10¹⁴ cm⁻², consistent with our previous observations (Indriolo et al. 2013b). The R(0) line shows a different profile from the R(1) and R(2) lines, suggesting that the excited states may be tracing a warmer gas component than the ground state. While less prominent than the HF features, absorption lines due to the P(34)–P(37) transitions of the v = 2–0 band of CO are also detected (Figure 1). These features are too weak to perform meaningful fits, but their depths are in agreement with predictions based on the temperature and column density of the warm CO component reported by Goto et al. (2019).

Zoom In Zoom Out Reset image size

Figure 8 shows portions of the spectra of AFGL 2136 IRS 1 and AFGL 4176 (A. Karska et al. 2020, in preparation), where the synthetic H₂O spectrum fails to reproduce five different absorption features. To determine the rest wavelengths of these features, we matched the line profiles for the unidentified lines with that of the HF R(0) line in the case of AFGL 4176, and that of the ${{\rm{H}}}_{3}^{+}$ $R{(1,1)}^{u}$ line (from Goto et al. 2019) in the case of AFGL 2136 IRS 1. The resulting rest wavelengths are 2.490070 μm, 2.490800 μm, 2.491885 μm, 2.492555 μm, and 2.493645 μm, with uncertainties of about 5 × 10⁻⁶ μm. We were unable to identify the one or more species that causes these absorption features, although based on the line profiles, the carrier must reside in the cooler foreground component, and not the warm gas that gives rise to H₂O absorption.

Zoom In Zoom Out Reset image size

Example line profiles from all detected species are presented in Figure 9. A variety of strong H₂O absorption lines are shown in panel (b), while weak H₂O lines are in the top half of panel (a). Panel (c) shows absorption out of the H₂O ground rotational state. Panel (d) shows absorption due to rovibrational transitions of ${{\rm{H}}}_{2}^{18}$O, NH₃, HCN, and C₂H₂, all of which have line profiles that are consistent with those observed for H₂O in the mid-infrared. The bottom half of panel (a) shows absorption due to HCl (Goto et al. 2013), CO, ${{\rm{H}}}_{3}^{+}$ (Goto et al. 2019), and HF, as well as two unidentified lines. The HCl P(4) absorption profile matches that of H₂O, while the profiles of ${{\rm{H}}}_{3}^{+}$ and the unidentified lines do not. HF is more complicated, as the R(0) line shows one component at the systemic velocity, while the R(1) and R(2) line profiles resemble the H₂O absorption. This suggests that the excited HF is tracing the same material as the H₂O, while the ground state HF is tracing a cooler foreground component that has been observed at T ≈ 60 K in the low-J transitions of the v = 2–0 band of CO (Goto et al. 2019). The NH₃, HCN, and C₂H₂ line profiles also suggest that these species reside in the warm gas component that gives rise to H₂O absorption.¹⁰ For comparison with the absorption lines, panel (e) shows the H₂O ${\nu }_{2}=1-1\,{5}_{\mathrm{5,0}}$–${6}_{\mathrm{4,3}}$ line at 232.687 GHz observed in emission by ALMA (extracted from the publicly available data cube provided by Maud et al. 2019).

Zoom In Zoom Out Reset image size

4.3.2. Emission Lines

In addition to the plethora of absorption features seen in the AFGL 2136 IRS 1 spectra, a few broad emission features are also present. These can be seen in Figure 1, and are due to the n = 19–5 and 18–5 transitions of the Pfund series of atomic H. The 17–5 transition was also seen in emission but suffered heavy interference from telluric H₂O lines and so was removed from the spectrum during the normalization process. A zoom in on the 19–5 and 18–5 emission features is shown in Figure 10. The 19–5 line is centered at 31 km s⁻¹ with a FWHM of 96 km s⁻¹, while the 18–5 line is centered at 39 km s⁻¹ with a FWHM of 127 km s⁻¹. Note, however, that our normalization and baseline removal procedures were optimized for the analysis of narrow absorption features, not broad emission features, so these values are highly uncertain.

Zoom In Zoom Out Reset image size

The CRIRES spectrum presented in Figure 1 is the average spectrum from observations taken on six different nights (see Table 1). Some unblended H₂O absorption lines were covered at all six epochs, so we can check for variability in the absorption profiles. Figure 11 shows the ν₃ 12_4,8–11_4,7 (top panel) and ν₃ 11_4,7–10_4,6 (bottom panel) transitions of H₂O from each individual night of observations. In the top panel, all of the observations from 2013 (39 days maximum separation) agree within the noise level, but the observation from 2012 suggests a deeper absorption feature with a less well-defined shoulder component toward longer wavelengths. In the bottom panel, all of the observations, including that from 2012, agree with each other within the noise level. However, this line in 2012 again lacks the well-defined shoulder that is present in all of the 2013 observations. It is difficult to draw any firm conclusions from these spectra. While the TEXES and EXES observations provide additional epochs in 2014 and 2016, respectively, due to the different instruments utilized, aperture sizes, and transitions observed, any such comparison would be highly uncertain. Still, the general agreement of the H₂O line profiles from the different telescopes in Figure 9 is indicative of no significant changes over a 4 yr timescale.

Zoom In Zoom Out Reset image size

The first indication of hot molecular gas near AFGL 2136 IRS 1 came from observations of the v = 1–0 band of ¹³CO and v = 2–1 band of ¹²CO in absorption (Mitchell et al. 1990). A hot component was determined to have $N{(}^{13}\mathrm{CO})=(2.5\pm 0.7)\times {10}^{17}$ cm⁻² and $T={580}_{-50}^{+60}$ K, containing about twice as much material as a cold component, and was hypothesized to arise in gas close to the central protostar. Recent observations of the v = 2–0 band of ¹²CO confirm this picture, finding $N{(}^{12}\mathrm{CO})=(2.8\pm 0.4)\times {10}^{19}$ cm⁻² and T = 534 ± 80 K (Goto et al. 2019). The absorption line profiles for ¹²CO v = 2–0 transitions out of highly excited levels (e.g., R(21) and R(23)) are consistent with those we find for H₂O, suggesting that both species are tracing the same component.

Observations with the Infrared Space Observatory Short Wavelength Spectrometer (ISO–SWS) have also revealed absorption from a variety of gas-phase molecules tracing warm gas, including CO₂, SO₂, H₂O, HCN, and C₂H₂. While the large wavelength coverage of ISO–SWS enables the observation and analysis of entire rovibrational bands, the resolution does not allow for the study of kinematics, individual line profiles, and state-specific column densities. Additionally, the determination of total column densities and excitation temperatures requires the use of assumed line widths. Analyses of molecular absorption bands observed toward AFGL 2136 IRS 1 by ISO–SWS have resulted in a range of excitation temperatures: ${T}_{\mathrm{ex}}({\mathrm{CO}}_{2})=300\pm 100$ K (Boonman et al. 2003), ${T}_{\mathrm{ex}}({\mathrm{SO}}_{2})={350}_{-50}^{+100}$ K (Keane et al. 2001a), ${T}_{\mathrm{ex}}({{\rm{H}}}_{2}{\rm{O}})={500}_{-150}^{+250}$ K (Boonman & van Dishoeck 2003), ${T}_{\mathrm{ex}}(\mathrm{HCN})={600}_{-50}^{+75}$ K (Lahuis & van Dishoeck 2000), and ${T}_{\mathrm{ex}}({{\rm{C}}}_{2}{{\rm{H}}}_{2})={800}_{-100}^{+150}$ K (Lahuis & van Dishoeck 2000). The cause of the large scatter in excitation temperatures is unclear, and it is difficult to evaluate if this scatter is real or if the excitation temperatures agree within uncertainties. Still, it is reasonable to associate these absorption features with the same (or a similar) gas component as that probed by the hot CO. Without kinematic information, though, it was not possible to confirm this hypothesis from the ISO–SWS observations alone.

In the past several years, observations of AFGL 2136 IRS 1 in the near- to mid-infrared have been made at much higher spectral resolution (∼3 km s⁻¹). Our previous observations of H₂O with CRIRES confirmed the high excitation temperature (${T}_{\mathrm{ex}}({{\rm{H}}}_{2}{\rm{O}})=506\pm 25$ K), but indicated a total H₂O column density about seven times larger than that determined from ISO–SWS ($N({{\rm{H}}}_{2}{\rm{O}})=(1.02\pm 0.02)\times {10}^{19}$ cm⁻² versus $N({{\rm{H}}}_{2}{\rm{O}})=(1.5\pm 0.6)\times {10}^{18}$ cm⁻²; Indriolo et al. 2013a; Boonman & van Dishoeck 2003). The resolved H₂O line profiles were in reasonable agreement with the hot CO component from Mitchell et al. (1990), but without access to the CO data, it was still not possible to directly compare the two species. Now, with our new H₂O observations and the recent CO v = 2–0 observations presented by Goto et al. (2019), we can definitively say that both species show the same absorption line profiles and very likely reside in the same physical component. This same component may also give rise to the NH₃, C₂H₂, and HCN absorption—as suggested by Figure 9—and this topic will be further addressed by A. G. Barr et al. (2020, in preparation). Although HCl shows a similar absorption profile to H₂O and CO, the excitation temperature inferred from the observed transitions is significantly lower: ${T}_{\mathrm{ex}}(\mathrm{HCl})\approx 250$ K (Goto et al. 2013). This is reminiscent of the scatter in the ISO–SWS results (if real), and potentially suggests that either (a) the species do not trace exactly the same material, but perhaps different layers of the same kinematic structure; or (b) excitation temperatures derived from individual molecules are not good indicators of the kinetic temperature of the gas since level populations may be influenced by radiative processes.

ALMA observations of AFGL 2136 IRS 1 have recently revealed emission from the hot molecular gas component as well. The H₂O ν₂ = 1–1, 5_5,0–6_4,3 transition at 232.687 GHz, which originates from a state E_u = 3461.9 K above ground, is observed tracing a compact disk (∼93 mas diameter) in Keplerian rotation about a central object with mass 45 ± 10 M_☉ (Maud et al. 2019). The velocity centroid along the disk minor axis is in agreement with the line-center velocity of the stronger component found from our Gaussian fitting procedure (24.8 km s⁻¹), suggesting that the H₂O absorption we observe in the infrared may arise in the disk now imaged at sub-mm wavelengths. Compact SiO emission (J = 5–4; 217.105 GHz; E_u = 31.3 K) has also been observed toward AFGL 2136 IRS 1 with ALMA. A portion of the SiO emission shows similar kinematics as H₂O and likely traces the disk (Maud et al. 2019), but emission near the systemic velocity is more spatially extended and hypothesized to arise in a rotating disk wind (Maud et al. 2018). Perhaps different conditions in these two components are responsible for the different rotation temperatures inferred from various molecules.

Even more compact than the disk around AFGL 2136 IRS 1 is a region producing emission in transitions indicative of ionized gas. The H30α recombination line is spatially unresolved with ALMA's 20 mas beam, but there is a hint of rotation in the same sense as the H₂O emission (Maud et al. 2019). It has a very broad line profile that is centered at 22.1 km s⁻¹, with FWHM ∼ 81.9 ± 1.7 km s⁻¹ (Maud et al. 2018). This is likely the same region that gives rise to the H Pfund n = 19–5 and 18–5 emission lines shown in our Figure 10, as the electrons cascade down toward the ground electronic state following recombination. Similar emission features from the n = 21–6 and 19–6 transitions of the Humphreys series of atomic H near 3.6 μm have also been observed toward AFGL 2136 IRS 1 (M. Goto 2020, private communication), and the hydrogen Brγ (n = 7–4) line has been observed with FWHM_Brγ ≈ 133 km s⁻¹ (Murakawa et al. 2013). Hydrogen Brγ emission in massive protostars has previously been interpreted as emission from an outflowing stellar wind or disk wind (e.g., Malbet et al. 2007; Murakawa et al. 2013), although the spatially unresolved H30α emission requires the emitting region to be much smaller than the circumstellar disk.

The optimal use of spectroscopy in studying massive protostars entails understanding what underlying characteristics cause the observed absorption and emission features. While the use of a simple isothermal absorbing slab model that is removed from and fully covering a background source can be instructive, its assumptions and limitations—many of which are described in Lacy (2013)—prevent such a model from reproducing observed spectra. To reiterate some of those points, (1) the central protostar will produce a temperature gradient in the surrounding gas and dust, (2) the absorbing gas can be mixed with the warm dust that serves as the "background" continuum, and (3) molecules both emit and absorb light. Here, we discuss how these and other effects can influence the observed spectra.

A covering fraction of less than one (i.e., clumpy absorbing material that does not cover the entire background source) is a scenario that has previously been invoked to explain observed infrared absorption spectra (e.g., Knez 2006; Knez et al. 2009; Barentine & Lacy 2012; Indriolo et al. 2015). In this case, saturated absorption lines will have a depth equal to the covering fraction, and all lines will appear weaker than expected, since total absorption corresponds to a relative intensity larger than zero. While some absorption features are about four times weaker than expected assuming the H₂O column density and excitation temperature inferred from CRIRES data alone (e.g., ν₂ 2_1,1–3_2,2 at 6.723374 μm), others are about as strong as expected (e.g., ν₂ = 2 − 1 7_1,6-7_2,5 at 6.726109 μm). Since covering fraction would affect all line depths equally, this explanation does not seem feasible. Even if the absorbing material does not fully cover the background source, the deepest absorption feature in the spectrum (H₂O ν₂ 3_2,1-3_1,2 at 6.075447 μm) reaches a transmission level of about 0.3, so a minimum covering fraction of 70% is required for the EXES observations.

One of the obvious differences among our spectra is the instruments used in performing the observations. In particular, the CRIRES, EXES, and TEXES observations were made using slits with different widths and lengths (widths of 02, 184, and 054, respectively), so light from slightly different regions is being observed in each case. The larger apertures could probe regions with different conditions than probed by the narrowest (CRIRES) slit from which we infer column density and temperature. Extended continuum emission (assuming the absorbing gas is not extended) or extended line emission could fill in the absorption features. The extended continuum explanation is equivalent to a covering fraction argument though, and so is ruled out. Extended line emission would have to be limited to only those transitions that are weaker than expected, but should also manifest as emission lines if a spectrum is extracted slightly off of the continuum peak. No H₂O line emission is seen when we extract a spectrum slightly off source, so this scenario seems unlikely as well.

Previous observations of molecular absorption and emission lines in the mid-infrared have demonstrated different behaviors of transitions within a rovibrational band. The ν₂ band of H₂O has shown some lines in emission and others in absorption (Gonzalez-Alfonso et al. 1998; González-Alfonso & Cernicharo 1999), as has the ν₂ band of NH₃ (Barentine & Lacy 2012). In both of these studies, though, there is a pattern that the R-branch lines appear in absorption, while the P-branch lines appear in emission, explained by the relative importance of the radiative pathways that populate and de-populate the vibrationally excited rotational levels. If this mechanism were responsible for the large scatter in the rotation diagrams—not to the point of causing emission, but just enough to fill in some of the absorption—then there should be a discrepancy in the column densities inferred from R-branch and P-branch transitions. Panel (c) of Figure 5 shows no correlation between "underpopulated" states and transition ΔJ though, so this emission mechanism is also unable to explain the observed spectra.

The parameter that shows the clearest correlation with a state being "underpopulated" is the depth of the absorption feature from which the state-specific column density was derived (Figure 6, panels (b), (d), and (f)). Excepting the absorption features observed with TEXES, deeper absorption lines tend to result in more "underpopulated" states. This can be explained by a scenario where the absorbing gas is well mixed with the warm dust that serves as a "background" source of continuum emission. Blackbody continuum emission predominantly comes from the τ ≈ 1 "surface," so any observed line absorption only comes from gas above that surface. The depth of the τ ≈ 1 "surface" depends on the dust opacity (and associated absorption coefficient, κ_ν), which is dependent on wavelength (Draine 2003) and on the distribution of dust grain sizes (e.g., Agurto-Gangas et al. 2019), so the layers of gas probed by observations at 2 μm, 6 μm, and 13 μm may be different. Because the absorption coefficient is also affected by line processes, κ_ν is larger at the wavelengths of transitions with larger spontaneous emission coefficients. This means that the physical depth where τ ≈ 1 is shallower for stronger lines than for weaker lines, such that there is a smaller column of gas in front of the "background" source. In this case, it is not that any levels are underpopulated with respect to LTE, just that they appear that way from a simple absorbing slab analysis since different transitions probe to different depths within the absorbing gas.

The provided scenario requires that the absorbing gas and source of continuum emission are well mixed. Mid-infrared (8–13 μm) interferometric observations of AFGL 2136 IRS 1 are suggestive of an elongated structure that has a size and orientation reasonably consistent with the recently reported sub-mm disk (de Wit et al. 2011; Boley et al. 2013; Maud et al. 2019). This indicates that the continuum emission at wavelengths longer than about 8 μm likely arises in the disk, and so it is reasonable to assume that the gas is mixed with the emitting dust for the TEXES observations. For silicate dust at the dust sublimation temperature (T ∼ 1300 K; Kobayashi et al. 2011), the blackbody spectrum peaks at about 2.2 μm, so the dust disk is expected to be bright in the near-infrared as well. Radiative transfer models of near-infrared interferometric observations of another massive protostar suggest that the dominant continuum emission source at 2.2 μm in that system is a circumstellar disk (Kraus et al. 2010, IRAS 13481−6124), so thermal emission from the warm dust disk may be the dominant continuum source in all of our observations from 2.5 to 13 μm (Beltrán & de Wit 2016).

Of the effects discussed previously, the mixing of warm dust and absorbing gas seems to be the most likely explanation for the observed H₂O absorption spectrum. ALMA observations have shown that the dust disk and gas disk have roughly the same physical extent (Maud et al. 2019), and it is plausible that the warm dust in the disk dominates the continuum at all observed wavelengths. The H₂O absorption line profiles, however, are not consistent with expectations in a scenario where the entire disk is responsible for continuum emission and molecular absorption. To better understand what may cause this discrepancy, a more detailed investigation of the line profiles is required.

Given the gas temperature inferred for the AFGL 2136 IRS 1 disk, line profiles are expected to be dominated entirely by orbital motions rather than thermal broadening. As such, information about gas kinematics is contained within the absorption and emission line profiles. H₂O absorption lines observed with CRIRES show two distinct absorption peaks while those observed with EXES and TEXES do not, but in all cases, lines are well described by two Gaussian components. The robustness of individual Gaussian fit parameters is difficult to assess, as in many cases intensity can be traded between the two components without adversely affecting the overall fit, so it is not particularly instructive to consider the properties of the two fit components separately. Indeed, the single temperature required to describe the rotation diagram for CRIRES data, and the relatively constant intensity ratio between the two components in the CRIRES data suggests that all of the absorption arises in a single physical component.

Absorption lines arising in circumstellar disks have previously been observed toward FU Ori type objects (young, low-mass stars experiencing luminosity bursts due to enhanced accretion; Hartmann & Kenyon 1996 and references therein). The fact that the lines appear in absorption rather than emission is indicative of a temperature inversion in the disk, where heating due to accretion causes the disk midplane to be warmer than the surface. Two properties observed in absorption lines toward FU Ori type objects are that (1) many absorption features are double-peaked; and (2) the width of absorption features tends to decrease with increasing wavelength. The double-peaked nature of the lines is explained by a disk in Keplerian rotation, with the two absorption peaks caused by the blueshifted and redshifted sides of the disk. While we do see double-peaked absorption lines toward AFGL 2136 IRS 1, the underlying explanation is not quite the same, as will be discussed later. The anticorrelation between line width and wavelength is explained by the temperature in the disk decreasing with increased distance from the central source. Continuum emission at shorter wavelengths preferentially arises in warmer material at small radii where orbital velocities are large, so the absorbing gas produces broad features. At longer wavelengths, the continuum emission arises from larger radii where orbital velocities are smaller and hence produce narrower absorption lines. An analysis of the distribution of full-width-at-zero-intensity (FWZI) line widths (defined as $| {v}_{2}-{v}_{1}| +3{\sigma }_{1}+3{\sigma }_{2}$) shows a marginal decrease (1.5 km s⁻¹ in the mean value) between H₂O lines at 2.5 μm and 5.7–6.7 μm (Figure 12). As FWZI line widths are broad (30 km s⁻¹ on average) and instrumental profiles have not been removed from the spectra, it is difficult to determine if this shift is real.

Zoom In Zoom Out Reset image size

The line profile for the 232 GHz H₂O emission extracted over the full extent of the disk is shown in Figure 9 panel (e). Inspection of the data cube channel by channel suggests that the flux at v ≲ 3 km s⁻¹ is an artifact of the image processing and that the H₂O emission from the disk extends from about 3 km s⁻¹ to about 45 km s⁻¹. Emission at these extreme velocities arises from the innermost portions of the disk along the major axis, and the velocity centroid along the disk minor axis is halfway between these two extremes at 24 km s⁻¹. The line profile has two emission peaks, with one each arising from the redshifted and blueshifted halves of the disk. All of these features are consistent with a disk in Keplerian rotation, as reported by Maud et al. (2019).

H₂O absorption profiles—and those of other molecules—clearly do not show the same characteristics as the 232 GHz H₂O emission. The strongest H₂O absorption occurs at about 25 km s⁻¹ with a second weaker component at about 33 km s⁻¹, and absorption extends from about 12 to 40 km s⁻¹. The narrower (with respect to emission), asymmetric line profiles seen in absorption suggest that if the H₂O absorption arises within the disk, it is not tracing the entire disk. In particular, the absence of absorption at the extreme velocities, and the lack of an absorption component at v ≤ 24 km s⁻¹, indicate that H₂O absorption probes neither the innermost portion of the disk along the major axis nor the blueshifted half of the disk. An examination of Figure 1(b) in Maud et al. (2019) provides a potential explanation for this observed behavior. There, the 1.3 mm continuum emission shows a clumpy structure within the disk after an assumed point source caused by free–free emission has been removed. The strongest continuum peaks within the disk are located along the minor axis, where v ≈ 24 km s⁻¹ is expected, and near the outer edge of the redshifted side, where v ≈ 33 km s⁻¹ is expected. These velocities are in good agreement with the two absorption components that are observed. If the clumpy structure seen in the 1.3 mm continuum emission is indicative of the 2–13 μm continuum structure, this could explain the molecular absorption profiles seen in the near- to mid-infrared.

As the ALMA observations were performed on 2017 October 4, the disk substructure at the time of the infrared observations would have been oriented slightly differently. Material at the inner and outer disk radii of 30 au and 120 au has orbital periods of about 24 yr and 196 yr, respectively (assuming a central mass of 45 M_☉), and would move by at most 18 mas and 9 mas in the plane of the sky over the 5.3 yr between our first CRIRES observation and the ALMA observations given a distance of 2.2 kpc. This suggests that the majority of the substructure identified in the ALMA observations was in a reasonably similar orientation at the time when the infrared observations were performed.

Finally, we note that there is no evidence for a high-velocity wind component in absorption toward AFGL 2136 IRS 1. CO absorption lines observed toward some massive protostars show broad blueshifted wings, extending over −100 km s⁻¹ from the systemic velocity (e.g., van der Tak et al. 1999). Neither the H₂O nor the CO (Mitchell et al. 1990; Smith et al. 2014; Goto et al. 2019) absorption lines observed toward AFGL 2136 IRS 1 show such a broad component. Most likely this is simply due to the orientation of the disk and outflow system, such that wind and outflow components are not probed by the pencil-beam line of sight toward the continuum source.

In our previous study of H₂O absorption observed with CRIRES toward AFGL 2136 IRS 1, we used a statistical equilibrium code to infer conditions within the absorbing gas (Indriolo et al. 2013a). The code modeled the excitation of the lowest 120 rotational states of ortho- and para-H₂O to determine the expected level populations under a variety of physical conditions. In this code, the effects of radiative trapping are treated using an escape probability method, and collisional rate coefficients between H₂O and H₂ were taken from Daniel et al. (2011) and were extrapolated to higher energy states for which calculations are unavailable using an artificial neural network method (Neufeld 2010). Radiative pumping effects in both rovibrational and rotational transitions due to the infrared continuum emission from AFGL 2136 IRS 1 were also included, assuming the spectral energy distribution implied by a fit to the data obtained by Murakawa et al. (2008)—particularly the TO model plotted in their figure 7. This analysis demonstrated that the observed state-specific column densities could be explained by purely collisional excitation, purely radiative excitation, or some contribution from both mechanisms between these two extremes. In the collisionally dominated case, the gas is at high density and temperature ($n({{\rm{H}}}_{2})\gtrsim 5\,\times \,{10}^{9}$ cm⁻³; T = 506 ± 25 K). In the radiatively dominated case, a T = 700 K blackbody subtends 2π sr (i.e., half of the sky) as seen by the gas, the gas temperature is unconstrained, and the density must be $n({{\rm{H}}}_{2})\lesssim {10}^{6}$ cm⁻³. Note that these two scenarios are limiting cases, and that some combination of collisional and radiative excitation is also viable, in which case the inferred constraints are relaxed (e.g., if excitation is primarily through collisions, but also due in some part to radiation, then the lower limit on the density will decrease compared to that in the purely collisional case).

We repeated this same statistical equilibrium analysis on the H₂O column densities presented here in Table 2 but were unable to find a model capable of simultaneously reproducing all values. This failure can be attributed to the large scatter in the rotation diagrams caused by different transitions probing different depths within the disk, as discussed previously. Limiting the analysis to only the CRIRES data produces the same results as in Indriolo et al. (2013a) that are summarized previously, so we do not present the entire analysis anew in this paper.

Given the limits on gas density and the measured H₂O column density, we can place constraints on the size of the region that contains H₂O. The relative abundance of H₂O with respect to H₂ in warm gas is expected to be about 10⁻⁴ (van Dishoeck et al. 2013 and references therein), so the radiatively dominated case suggests an upper limit of $n({{\rm{H}}}_{2}{\rm{O}})\lesssim 100$ cm⁻³. With $N({{\rm{H}}}_{2}{\rm{O}})=8.25\,\times \,{10}^{18}$ cm⁻², this requires a minimum path length of l ≳ 8.25 × 10¹⁶ cm or about 5500 au. The collisionally dominated case requires $n({{\rm{H}}}_{2}{\rm{O}})\gtrsim 5\,\times \,{10}^{5}$ cm⁻³, corresponding to a maximum path length of $l\lesssim 1.65\times {10}^{13}$ cm, or about 1.1 au. It is clear that the path length over which absorption occurs in the radiatively dominated case is incompatible with the r ∼ 120 AU size scales of the AFGL 2136 IRS 1 disk. The path length in the collisionally dominated case, however, is consistent with the picture where the observed H₂O absorption is coming from the upper layers of the disk. That said, the fact that the H₂O resides within the disk means that blackbody radiation from the dust disk must have some influence on the excitation of H₂O, so the absorbing layer can be thicker than that determined for the scenario where only collisions are considered.

Although Figure 11 shows the possibility of marginal variability in H₂O absorption line profiles, the similarity of the line profiles in the CRIRES, EXES, and TEXES data shown in Figure 9 suggests no significant changes in the absorbing material over 4 yr timescales. This is further evidence against a scenario where the background continuum source and absorbing gas are separated (e.g., gas in the disk absorbing light from the central protostar), because the absorbing material in a pencil beam toward the central object will change as it orbits. The time between our first CRIRES observation and last EXES observation was 1357 days. Over this time period, gas at radii of 30 au and 120 au in Keplerian rotation about a 45 M_☉ star would move 28.6 au and 14.3 au along their respective orbits. The disk material would have to be uniform over these size scales in order for the absorption profile to remain nearly constant in a pencil beam toward the central object over 4 yr. If portions of the disk itself serve as the background continuum source though, then the material above the τ ≈ 1 "surface" will be probed by absorption spectroscopy, regardless of where the material is in its orbit. Time variability in line profiles is still expected in this case though, as the 1.3 mm continuum emission shows substructure (Maud et al. 2019). As regions of brighter continuum emission orbit the central object, the absorption from the gas associated with these regions will shift in velocity.

In our observations, absorption due to both the ${{\rm{H}}}_{2}^{16}{\rm{O}}$ and ${{\rm{H}}}_{2}^{18}{\rm{O}}$ isotopologues of H₂O is detected. With only two unblended ${{\rm{H}}}_{2}^{18}{\rm{O}}$ lines observed, we cannot estimate the total column density of that isotopologue, but we can compare state-specific column densities. Transitions probing the ${J}_{{K}_{a},{K}_{c}}={3}_{\mathrm{1,2}}$ state of both ${{\rm{H}}}_{2}^{16}{\rm{O}}$ and ${{\rm{H}}}_{2}^{18}{\rm{O}}$ are observed in our spectrum, providing column densities of $N({{\rm{H}}}_{2}^{16}{\rm{O}}| {3}_{\mathrm{1,2}})=(21.1\pm 0.84)\times {10}^{15}$ cm⁻² and $N({{\rm{H}}}_{2}^{18}{\rm{O}}| {3}_{\mathrm{1,2}})=(2.23\pm 0.86)\times {10}^{15}$ cm⁻² (see Table 2). This state-specific analysis results in ${{\rm{H}}}_{2}^{16}{\rm{O}}/{{\rm{H}}}_{2}^{18}{\rm{O}}=9.5$, much smaller than the oxygen isotopic ratio found in the solar system (¹⁶O/¹⁸O ≈ 500; Asplund et al. 2009 and references therein). However, the ${{\rm{H}}}_{2}^{16}{\rm{O}}$ ${\nu }_{2}=1-0\,{3}_{\mathrm{2,1}}$–3_1,2 transition is one of the strongest absorption features in our spectrum, suggesting that it traces a much smaller region than the optically thin ${{\rm{H}}}_{2}^{18}{\rm{O}}$ ${\nu }_{2}=1-0\,{3}_{\mathrm{2,1}}$–3_1,2 transition. Rather than representing the real isotopic ratio then, this is another indicator of the optical depth effects discussed previously. The best-fit temperature and ${{\rm{H}}}_{2}^{16}{\rm{O}}$ column density can be used to estimate the state-specific column densities for ${{\rm{H}}}_{2}^{16}{\rm{O}}$ corresponding to the two observed ${{\rm{H}}}_{2}^{18}{\rm{O}}$ transitions for use in determining isotopic ratios. Using this method, we find $N({{\rm{H}}}_{2}^{16}{\rm{O}}| {1}_{\mathrm{0,1}})=1.79\,\times \,{10}^{17}$ cm⁻² and $N({{\rm{H}}}_{2}^{16}{\rm{O}}| {3}_{\mathrm{1,2}})=2.71\,\times \,{10}^{17}$ cm⁻². Together with the measured ${{\rm{H}}}_{2}^{18}{\rm{O}}$ state-specific column densities, these values give ${{\rm{H}}}_{2}^{16}{\rm{O}}/{{\rm{H}}}_{2}^{18}{\rm{O}}\approx 206$ and ${{\rm{H}}}_{2}^{16}{\rm{O}}/{{\rm{H}}}_{2}^{18}{\rm{O}}\approx 121$. Although still below the solar system isotopic ratio, these are well within the wide range of values found in other stellar atmospheres (Hinkle et al. 2016).