# MINDS: The DR Tau disk

## II. Probing the hot and cold H₂O reservoirs in the JWST-MIRI spectrum

Milou Temmink¹, Ewine F. van Dishoeck¹,², Danny Gasman³, Sierra L. Grant², Benoît Tabone⁴,
Manuel Güdel⁵,⁶, Thomas Henning⁷, David Barrado⁸, Alessio Caratti o Garatti⁹,¹⁰, Adrian M. Glauser⁶,
Inga Kamp¹¹, Aditya M. Arabhavi¹¹, Hyerin Jang¹², Nicolas Kurtovic², Giulia Perotti⁷,
Kamber Schwarz⁷, and Marissa Vlasblom¹

*(Affiliations can be found after the references)*

Received 12 April 2024 / Accepted 5 July 2024

## ABSTRACT

**Context.** The Medium Resolution Spectrometer (MRS) of the Mid-InfraRed Instrument (MIRI) on the James Webb Space Telescope (JWST) gives insights into the chemical richness and complexity of the inner regions of planet-forming disks. Several disks that are compact in the millimetre dust emission have been found by Spitzer to be particularly bright in H₂O, which is thought to be caused by the inward drift of icy pebbles. Here, we analyse the H₂O-rich spectrum of the compact disk DR Tau using high-quality JWST-MIRI observations.

**Aims.** We infer the H₂O column densities (in cm⁻²) using methods presented in previous works, as well as introducing a new method to fully characterise the pure rotational spectrum. We aim to further characterise the abundances of H₂O in the inner regions of this disk and its abundance relative to CO. We also search for emission of other molecular species, such as CH₄, NH₃, CS, H₂, SO₂, and larger hydrocarbons; commonly detected species, such as CO, CO₂, HCN, and C₂H₂, have been investigated in our previous paper.

**Methods.** We first use 0D local thermodynamic equilibrium (LTE) slab models to investigate the excitation properties observed in different wavelength regions across the entire spectrum, probing both the ro-vibrational and rotational transitions. To further analyse the pure rotational spectrum (≥10 µm), we use the spectrum of a large, structured disk (CI Tau) as a template to search for differences with our compact disk. Finally, we fit multiple components to characterise the radial (and vertical) temperature gradient(s) present in the spectrum of DR Tau.

**Results.** The 0D slab models indicate a radial gradient in the disk, as the excitation temperature (emitting radius) decreases (increases) with increasing wavelength, which is confirmed by the analysis involving the large disk template. To explain the derived emitting radii, we need a larger inclination for the inner disk (i ~ 10–23°), agreeing with our previous analysis on CO. From our multi-component fit, we find that at least three temperature components (T₁ ~800 K, T₂ ~470 K, and T₃ ~180 K) are required to reproduce the observed rotational spectrum of H₂O arising from the inner R<sub>em</sub> ~0.3–8 au. By comparing line ratios, we derived an upper limit on the column densities (in cm⁻²) for the first two components of log₁₀(N) ≤18.4 within ~1.2 au. We note that the models with a pure temperature gradient provide as robust results as the more complex models, which include spatial line shielding. No robust detection of the isotopologue H₂ ¹⁸O can be made and upper limits are provided for other molecular species.

**Conclusions.** Our analysis confirms the presence of a pure radial temperature gradient present in the inner disk of DR Tau, which can be described by at least three components. This gradient scales roughly as ~R<sub>em</sub><sup>-0.5</sup> in the emitting layers, in the inner 2 au. As the observed H₂O is mainly optically thick, a lower limit on the abundance ratio of H₂O/CO~0.17 is derived, suggesting a potential depletion of H₂O. Similarly to previous work, we detect a cold H₂O component (T ~ 180 K) originating from near the snowline, now with a multi-component analysis. Yet, we cannot conclude whether an enhancement of the H₂O reservoir is observed following radial drift. A consistent analysis of a larger sample is necessary to study the importance of drift in enhancing the H₂O abundances.

**Key words.** astrochemistry – protoplanetary disks – stars: variables: T Tauri, Herbig Ae/Be – infrared: general

## 1. Introduction

H₂O is a key ingredient in making habitable planets. As the bulk of exoplanets form in the inner-most regions (<10 au; Morbidelli et al. 2012; Dawson & Johnson 2018) of planet-forming disks, their (atmospheric) composition is determined by the available elemental abundances. Therefore, it is of importance to characterise these regions in detail and, in particular, to analyse the H₂O emission present.

The recently launched James Webb Space Telescope (JWST; Rigby et al. 2023) and, in particular, the Medium Resolution Spectroscopy mode (MRS; Wells et al. 2015) of the Mid-Infrared Instrument (MIRI; Rieke et al. 2015; Wright et al. 2015), provides the best opportunity to study the chemical composition

of these inner-most regions of planet-forming disks. The wide wavelength range of JWST-MIRI (4.9–27.9 µm) covers many transitions of H₂O. These include the ro-vibrational lines at the shorter wavelengths between ~5.0 and 8.0 µm and the rotational transitions longward of 10.0 µm (Meijerink et al. 2009). This forest of H₂O lines across the mid-infrared wavelength provides insights into the different parts of the inner disk: the lines at shorter wavelengths are thought to probe the innermost gas, whereas those at longer wavelengths probe colder gas located at larger radial distances (e.g. Blevins et al. 2016; Banzatti et al. 2017, 2023b, and Gasman et al. 2023). Observations suggest that H₂O vapour is prevalent in the spectra of T-Tauri disks, even in those with inner, dust-depleted cavities, albeit with lower line strengths (Perotti et al. 2023; Schwarz et al. 2024).
Furthermore, Banzatti et al. (2020) have found a correlation between the flux of strong H₂O lines measured with the *Spitzer* Space Telescope and the radial dust disk size ($R_{\text{dust}}$) as measured with ALMA. The proposed explanation for the correlation is driven by the H₂O abundances expected for three types of disks: compact disks ($R_{\text{dust}} \le 60$ au) with efficient radial drift; large disks ($60 \le R_{\text{dust}} \le 300$ au) with substructures; and large disks with substructures and an inner cavity. Substructures are able to trap icy pebbles in the outer disk, preventing them from drifting inside the H₂O snowline and enhancing the H₂O column density in the inner disk. Subsequently, these large disks are thought to have low H₂O column densities in the inner disk, which are expected to be further depleted in the presence of an inner cavity. For the compact dust disks, where radial drift is thought to be very efficient (e.g., Facchini et al. 2019) and substructures are found to be less common (Long et al. 2019), the H₂O column density is thought to be enhanced through the sublimation of H₂O-ice. However, not all compact disks show strong H₂O emission, such as DN Tau, for example (Pontoppidan et al. 2010; Salyk et al. 2011; Banzatti et al. 2020).

In this work, we analyse the H₂O emission in the JWST-MIRI spectrum of DR Tau, a compact T-Tauri (K6) disk (<60 au, Long et al. 2019) located at a distance of ~195 pc (Gaia Collaboration 2018) in the Taurus star-forming region. DR Tau has a mass of M = 0.93 M⊙, an effective temperature of T<sub>eff</sub> = 4205 K, and a luminosity of L = 0.63 L⊙ (Long et al. 2019). Observations with *Spitzer* have shown that DR Tau has one of the highest line-to-continuum ratios and contains a rich H₂O reservoir (e.g., Salyk et al. 2011 and Banzatti et al. 2020). Using ground-based observatories, various previous studies have analysed some of the bright H₂O lines at much higher spectral resolution at both near- and mid-infrared wavelengths (Salyk et al. 2008; Najita et al. 2018; Salyk et al. 2019; Banzatti et al. 2023b). High-resolution, ground-based ro-vibrational CO observations at ~4.6–5.3 µm are also available for this disk, which have been analysed in Bast et al. (2011); Brown et al. (2013); Banzatti et al. (2022) and Temmink et al. (2024).

Temmink et al. (2024) used rotational diagrams to investigate the excitation properties of the various CO isotopologues: ¹²CO, ¹³CO, and C¹⁸O. For the optically thin isotopologue C¹⁸O an excitation temperature of T ~ 975 K, a column density of N ~2.0×10¹⁶ cm⁻², and an emitting radius of R<sub>CO</sub> ~ 0.23 au were derived. Using these parameters and accounting for the involved isotopologue ratio, a total number of molecules of N<sub>CO</sub> ~ 4.1 × 10⁴⁴ was found for CO. The physical parameters derived from the high spectral resolution CO observations were shown to also provide an explanation for the CO emission observed with JWST-MIRI. In addition, Temmink et al. (2024) investigated the emission of CO₂, HCN, and C₂H₂. Left unanalysed, however, was the plethora of H₂O lines present in the spectrum. Due to the strength of the H₂O lines, the high line-to-continuum ratio, and the disk having the most ground-based, ro-vibrational and rotational H₂O transitions observed to date (see, for example, Najita et al. 2018; Salyk et al. 2019; Banzatti et al. 2023b), DR Tau is one of the best candidates for an extensive analysis of the H₂O emission and to make a detailed comparison with the physical structure and abundance of CO, under the assumption that these molecules are at least partially co-spatial. In addition, due to its compactness, DR Tau is also a good candidate to study the effects of radial drift on the observable H₂O abundance.

In this paper, we provide a deep analysis of the H₂O-rich spectrum using three different methods: we analyse various wavelength regions in the spectrum using 0D local thermal

equilibrium (LTE) slab models. Additionally, we investigate the difference between the spectrum of a compact disk with that of a large, structured disk. Finally, we introduce a multi-component, radial and vertical gradient slab model fitting technique, effectively providing a 1D LTE slab model, to derive the temperature and column density gradients of the emitting layers in the inner disk.

This paper is structured as follows: in Sect. 2 we describe the JWST-MIRI observations of DR Tau. Section 3 provides the used analysis methods and the accompanying results, which are further discussed in Sect. 4. Finally, our conclusions are summarised in Sect. 5.

## 2. Observations

DR Tau has been observed with JWST-MIRI as part of the JWST Guaranteed Time Observations Program MIRI mid-INfrared Disk Survey (MINDS, PID: 1282, PI: T. Henning; Henning et al. 2024). The details of the observations are described in Temmink et al. (2024). For the analysis in this work, we have re-reduced the observations using a standard pipeline reduction (Version 1.12.3, Bushouse et al. 2023) and the photom updates released in November 2023, which significantly improve the reduction of Channel 4. This reduction uses the residual fringe correction implemented in the pipeline, uses an annulus background subtraction, and the spectrum is extracted using an aperture with a size of 2x the full width at half maximum (FWHM). We note that the photom updates do not have a significant impact on the other channels, indicating that the results listed in Temmink et al. (2024) still hold.

We have used the continuum subtraction method as described in Temmink et al. (2024), which includes Channel 4 now as well. This continuum subtraction method first filters out downwards spikes before estimating the baseline with the 'Iterative Reweighted Spline Quantile Regression' (IRSOR) method included in the python-package PYBASELINES (Erb 2022). The estimated continuum is also visible in Fig. 1.

Our results, the multi-component slab model fits (see Sect. 3.3) and the potential (non-)detection of H₂ ¹⁸O (see Sect. 3.5) may be partially influenced by the wavelength calibration of JWST-MIRI/MRS. The wavelength accuracy, as shown by Argyriou et al. (2023) and Pontoppidan et al. (2024), can be off by 40–90 km s⁻¹ and the offset is the largest in Channel 4 as clearly shown in Fig. 5 of Pontoppidan et al. (2024). As Channel 4 covers a large range of important H₂O transitions, as well as the strongest H₂ ¹⁸O transitions observable with JWST-MIRI, the velocity shifts can significantly impact the slab models and/or one's ability of detecting H₂ ¹⁸O.

To account for the wavelength correction, we have used, similar to Pontoppidan et al. (2024), a plethora of CO and H₂O lines to estimate the offset between the spectrum and an initial model through a cross-correlation technique. In our estimates, we have accounted for the heliocentric velocity of DR Tau ($v_{\text{hel}}$ ~ 27.64 km s⁻¹; Ardila et al. 2002). In total, we have used 483 (blended) transitions, yielding corrections between ~−100 km s⁻¹ and ~120 km s⁻¹ with the largest offsets found in Channel 4. As the silicate feature is dominated by noise and no strong H₂O transitions are found between ~8 and ~10 µm, we have not determined the offsets for any possible transitions present in this region. In addition, as the offsets for the shortest wavelengths are mostly concentrated around zero, we have also not applied corrections for these regions. The corrections are thus only applied for the longer wavelengths (>10 µm), where
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Full JWST-MIRI spectrum (4.9–27.5 µm) of DR Tau. The different wavelength ranges (subbands) of each MIRI/MRS channel are shown in blue ('A'), green ('B'), and orange ('C'), respectively. The wavelength regions over which H₂O will be analysed are indicated by the horizontal bars. The continuum fit is indicated by the black line.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Estimated wavelength corrections (black circles) for the selected lines. The red lines are the polynomials of degree 3 fitted for each JWST-MIRI subband and are used in correcting the wavelength calibration.</figcaption></figure>

the offsets are found to be the largest. Figure 2 displays our offsets found for the different CO and H₂O transitions.

To properly determine the wavelength corrections across the spectrum of DR Tau, we have fitted polynomials of degree 3 through the found offsets for each JWST-MIRI subband. These polynomials, indicated by the red lines in Fig. 2, have been used to correct the spectrum. We note that the corrections of Pontoppidan et al. (2024) have been implemented for the longest wavelengths in a newer version of the JWST reduction pipeline than used in this work (priv. comm.).

**3. Analysis and results**

We analysed the H₂O emission in DR Tau using four different methods: first, in Sect. 3.1, similar to Gasman et al. (2023), we investigated the H₂O emission across separate regions in the spectrum using 0D local-thermal equilibrium (LTE) slab models (see Sect. 3.1). The slab model fitting procedure is described in detail in Temmink et al. (2024) (see also Grant et al. 2023 and Tabone et al. 2023). As the new reduction yields slightly different fluxes, especially in Channel 4, we report in Table A.1 new values for the median flux, the signal-to-noise ratio (S/N) as given by the JWST Exposure Time Calculator (ETC)<fnref n="1" />, which

accounts for the degrading transmission<fnref n="2" /> of JWST-MIRI/MRS, and for the estimated noise on the continuum used in our fitting procedure. For the H₂O slab models, we used grids where the column density (N in cm⁻<fnref n="2" />) was allowed to vary between 12 and 22 in log₁₀-spacing using a of ΔN = 0.1. We note that these column densities probe the upper layers of the disk, above the height where the dust becomes optically thick, or, in the case of optically thick lines, above the height where the emission lines become optically thick (Bruderer et al. 2015; Woitke et al. 2018). The temperature (T) was allowed to vary between 150 and 2500 K in steps of ΔT = 25 K, whereas the emitting radius (R<sub>em</sub>) was allowed to take on values between 0.01 and 10 au, using step sizes of ΔR = 0.024 au. The emitting radius is used as a parametrisation of the emitting area, assuming A = πR<sub>em</sub><fnref n="2" />. We note that the emitting area does not have to be a circle, it can have any shape (i.e. an annulus) at any radial distance from the star. In addition to regions listed in Gasman et al. (2023) (5.0–6.5, 13.6–16.3, 17.0–23.0, and 23.0–27.0 µm), we also investigate the H₂O emission in the regions between 7.0–9.5 and 10.5–13.0 µm. The different wavelength regions are also indicated in Fig. 1.

Second (see Sect. 3.2), we followed the method described in Banzatti et al. (2023a), where a continuum-subtracted spectrum of a large disk (CI Tau) was used to investigate a second, colder H₂O component in a compact disk (GK Tau), after scaling for the differences in distances and luminosity of the H₂O lines with highest upper energy levels (6000–8000 K). They expect this second component to trace the enrichment of the H₂O gas near the snowline, following the sublimation of icy mantles of drifting grains. As DR Tau is a compact disk, we expect pebble drift to be efficient and, subsequently, we expect such a second H₂O component should also be visible. Here, we also used the spectrum of CI Tau as the template spectrum to investigate the presence of an excess, cold second component in DR Tau (see Sect. 3.2.2). For consistency, the CI Tau observations have been reduced in the same way as done for DR Tau (see Sect. 2).

As a third method (Sect. 3.3), we investigate how well 0D LTE slab models can reproduce the observed spectrum when implementing radial and vertical temperature gradients.

<footnote n="1">¹ ETC: https://jwst.etc.stsci.edu/</footnote>

<footnote n="2">² For more information on the degrading transmission, see: https://www.stsci.edu/contents/news/jwst/2023/miri-mrs-reduced-count-rate-update</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Inferred gas temperature (left) and emitting radii (middle) for the different wavelength regions, and the emitting temperature as a function of the emitting radius (right). We show the results for the slab model fits without (squares) and with (circles) line overlap. The colours indicate the varying gas temperatures, whereas the horizontal bars denote the wavelength ranges of the regions. Errorbars are shown for the derived quantities, where those for the models without line overlap are shown in black and those for the models with line overlap are shown in red. In addition, we show for comparison the results (blacks stars) of the three component fit (approach 3, see Sect. 3.3.3) in the right panel.</figcaption></figure>

Finally, we obtained limits on the H₂O column density (in cm⁻²) by comparing the ratios of the Einstein-A coefficients and the line fluxes of H₂O line pairs with the same upper energy level (Sect. 3.4, see also Gasman et al. 2023). We also searched for emission of the rare isotopologue H₂ ¹⁸O in Sect. 3.5, using the isolated lines listed in Calahan et al. (2022).

Additionally, in Sect. 3.6, we fitted an LTE slab model to the OH emission over the wavelength region between 10.0 and 27.5 µm and investigated the presence of other molecular species, such as CH₄, NH₃, and larger hydrocarbons. We compare their emission properties with those of H₂O, CO₂, HCN, and C₂H₂, of which the final three are presented in Temmink et al. (2024).

### 3.1. Single H₂O LTE slab models

The JWST-MIRI observations of DR Tau show a plethora of both ro-vibrational (<10 µm) and rotational (>10 µm) transitions of H₂O. Following Gasman et al. (2023), we investigated the excitation properties of H₂O across different wavelength regions, covering both ro-vibrational and rotational transitions, using OD-LTE slab models. We have used slab models without and with line overlap, i.e. mutual shielding of H₂O lines (i.e. no shielding by other molecules, see Tabone et al. 2023 for a full description). Table 1 lists the best fitting slab model parameters for both the models without line overlap (top half) and with line overlap (bottom half). Figures displaying the best fitting slab models without and with line overlap are available online<fnref n="3" />.

The best fitting models were determined using a reduced χ²-minimisation method, following, for example, Grant et al. (2023), Tabone et al. (2023), and Temmink et al. (2024). The reduced χ² is determined using the following equation:

$$ \chi_{\text{red}}^2 = \frac{1}{N_{\text{obs}}} \sum_{i=1}^{N_{\text{obs}}} \frac{(F_{\text{obs},i} - F_{\text{mod},i})^2}{\sigma^2}. \quad (1) $$

Here, $N_{\text{obs}}$ is the number of spectral resolution elements, covering isolated H₂O transitions, used in the fitting, $F_{\text{obs}}$ and $F_{\text{mod}}$ are the corresponding observed and model flux, respectively, and $\sigma$ denotes the uncertainty on the flux (see Table A.1). The uncertainties on the best-fit parameters are taken from the confidence

<table><caption>**Table 1.** Best fit parameters of the H₂O slab models without (top) and with line overlap (bottom).</caption><thead><tr><th>Region (µm)</th><th>log<sub>10</sub> (N)<sup>(a)</sup></th><th>T (K)</th><th>R (au)</th><th>N</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">Without line overlap</td></tr><tr><td>5.0–6.5</td><td>18.4<sup>+0.2</sup><sub>-0.1</sub></td><td>1000<sup>+250</sup><sub>-50</sub></td><td>0.21<sup>+0.02</sup><sub>-0.04</sub></td><td>7.79×10<sup>43</sup></td></tr><tr><td>7.0–9.5</td><td>18.3<sup>+0.1</sup><sub>-0.1</sub></td><td>1000<sup>+50</sup><sub>-100</sub></td><td>0.21<sup>+0.02</sup><sub>-0.00</sub></td><td>6.19×10<sup>43</sup></td></tr><tr><td>10.5–13.0</td><td>18.1<sup>+0.3</sup><sub>-1.8</sub></td><td>850<sup>+50</sup><sub>-50</sub></td><td>0.87<sup>+5.18</sup><sub>-0.20</sub></td><td>6.70×10<sup>44</sup></td></tr><tr><td>13.6–16.3</td><td>18.4<sup>+0.0</sup><sub>-0.0</sub></td><td>675<sup>+0</sup><sub>-0</sub></td><td>0.93<sup>+0.00</sup><sub>-0.00</sub></td><td>1.53×10<sup>45</sup></td></tr><tr><td>17.0–23.0</td><td>18.9<sup>+0.6</sup><sub>-0.3</sub></td><td>450<sup>+25</sup><sub>-50</sub></td><td>1.31<sup>+0.14</sup><sub>-0.12</sub></td><td>9.58×10<sup>45</sup></td></tr><tr><td>23.0–27.0</td><td>20.1<sup>+2.0</sup><sub>-2.1</sub></td><td>375<sup>+350</sup><sub>-175</sub></td><td>1.55<sup>+1.92</sup><sub>-0.64</sub></td><td>2.13×10<sup>47</sup></td></tr><tr><td colspan="5" style="text-align: center;">With line overlap</td></tr><tr><td>5.0–6.5</td><td>18.5<sup>+0.1</sup><sub>-0.0</sub></td><td>1150<sup>+100</sup><sub>-75</sub></td><td>0.17<sup>+0.00</sup><sub>-0.00</sub></td><td>6.43×10<sup>43</sup></td></tr><tr><td>7.0–9.5</td><td>18.5<sup>+0.2</sup><sub>-0.1</sub></td><td>850<sup>+75</sup><sub>-75</sub></td><td>0.25<sup>+0.02</sup><sub>-0.02</sub></td><td>1.39×10<sup>44</sup></td></tr><tr><td>10.5–13.0</td><td>18.0<sup>+0.4</sup><sub>-1.3</sub></td><td>850<sup>+50</sup><sub>-50</sub></td><td>0.95<sup>+2.86</sup><sub>-0.30</sub></td><td>6.35×10<sup>44</sup></td></tr><tr><td>13.6–16.3</td><td>18.5<sup>+0.0</sup><sub>-0.0</sub></td><td>675<sup>+25</sup><sub>-0</sub></td><td>0.89<sup>+0.00</sup><sub>-0.06</sub></td><td>1.76×10<sup>45</sup></td></tr><tr><td>17.0–23.0</td><td>18.5<sup>+0.3</sup><sub>-0.2</sub></td><td>500<sup>+25</sup><sub>-50</sub></td><td>1.35<sup>+0.14</sup><sub>-0.08</sub></td><td>4.05×10<sup>45</sup></td></tr><tr><td>23.0–27.0</td><td>19.5<sup>+2.6</sup><sub>-1.8</sub></td><td>400<sup>+350</sup><sub>-200</sub></td><td>1.69<sup>+2.36</sup><sub>-0.68</sub></td><td>6.35×10<sup>46</sup></td></tr></tbody></table>

Notes. (a): N is given in units of in cm⁻².

intervals, which are determined for, respectively, 1σ, 2σ, and 3σ as $\chi^2_{\text{red,min}} + 2.3$, $\chi^2_{\text{red,min}} + 6.2$ and $\chi^2_{\text{red,min}} + 11.8$ (Avni 1976; Press et al. 1992).

In accordance with previous works (e.g. Blevins et al. 2016; Banzatti et al. 2023a; Gasman et al. 2023), the slab models probe the inner, hotter regions at the shortest wavelengths, whereas the outer, colder regions are traced by the longer wavelengths. The trends hold for both the models with and without line overlap. Figure 3 visualises the decreasing (increasing) trend of temperature (emitting radius) with wavelength for the models without overlap (squares) and those with line overlap (circles). We find that the temperatures decrease from 1000 K to 400 K over a span of 0.2 to 2.0 au, if $R_{\text{em}}$ corresponds to the actual radius of the emitting area. Additionally, we show the temperature as a function of radius in the right panel of Fig. 3. This panel also includes the results from our three component fit (using method 3, see Sect. 3.3) and the 1σ uncertainties.

All fits yield generally well constrained values for the column densities, excitation temperatures and emitting areas,

<footnote n="3">³ https://zenodo.org/records/12703820: Figs. B1-1 to B1-6 display the slab models without line overlap, whereas Figs. B2-1 to B2-6 display those with line overlap.</footnote>
parameterised by the emitting radius $R_{em}$. The uncertainties, which are mostly within ±0.5 for the logarithm of the column density (in cm⁻²), ±100 K for the excitation temperature, and ±0.20 au for the emitting radius, are the largest for the fits in region 6 (23.0–27.0 µm), likely the effect of this wavelength region having the lowest sensitivity. In addition, the upper uncertainty on the emitting radius for the third region is found to be higher compared to the other uncertainties. This is likely the result of this region having fewer bright, optically thick lines. Using the results for the inner two regions, those traced by the ro-vibrational H₂O transitions, we can make a comparison with those for CO as analysed in Temmink et al. (2024), since those arise from a similar region of ~0.2 au. This comparison is made in Sect. 4.4.

We find that our slab models to the JWST-MIRI spectrum of DR Tau yield very similar results for the excitation temperatures and emitting radii for both the vibrational and rotational transitions of H₂O compared to the previously mentioned ground-based, high spectral resolution observations (Najita et al. 2018; Salyk et al. 2019; Banzatti et al. 2023b). Furthermore, Banzatti et al. (2023b) provides linewidths for the vibrational H₂O lines at 5.0 µm (FWHM ~ 22 km s⁻¹) and the rotational transitions at 12.4 µm (FWHM ~ 16 km s⁻¹), so we can estimate the emitting radii from these line profiles using Kepler's law: $R = GM_* (\sin(i)/\Delta v)^2$, where i is the inclination of the disk (i ~ 5.4, as derived from ALMA observations; Long et al. 2019) and Δv is the linewidth, here taken equal to half the width at half maximum (HWHM). We derive emitting radii of $R_{em} \sim 0.06$ au and ~0.11 au for the 5.0 and 12.4 µm high-resolution transitions, respectively. Both of these are significantly smaller than the derived emitting radii from our JWST observations and the slab models fitted to the ground-based, high-resolution observations. A similar discrepancy was found in our previous work, analysing the high spectral resolution CO observations (Temmink et al. 2024). There we found that the emitting radii derived from the slab models and those using the linewidths for CO agree if the (inner) disk inclinations are increased to $i_{inner} \sim 10–23°$. Notably, observations of the VLTI-GRAVITY instrument suggest at similar inclinations for the inner disk ($i_{inner} \sim 18^{+10}_{-18}°$; GRAVITY Collaboration 2021). Using these inclinations, we derive emitting radii in the ranges of $R_{em} \sim 0.21–0.96$ au for vibrational transitions at 5 µm and $R_{em} \sim 0.39–1.81$ au, which do agree with the emitting radii derived from the slab models. As mentioned in Temmink et al. (2024), these derivations hold for the assumption that the emitting radii correspond to a circular region enclosing the host star and may suggest a misalignment between the inner and outer disk of DR Tau.

## 3.2. Cold and warm H₂O reservoirs: large disk template

### 3.2.1. Identifying unblended lines

Evidence for the presence of multiple rotational H₂O components in the JWST-MIRI spectra was introduced by Banzatti et al. (2023a) and Gasman et al. (2023). In a similar fashion to Banzatti et al. (2023a), we first identified the unblended, high-energy (6000 ≤ $E_{up}$ ≤ 8000 K) H₂O lines, which all happen to have Einstein-A coefficients $A_{ul} > 10$ s⁻¹. Figures displaying $A_{ul}$ as a function of $E_{up}$ (left panel) and wavelength (right panel) are available online<fnref n="4" />. As our analysis includes many more H₂O lines

compared to Banzatti et al. (2023a), we have identified fewer isolated lines. This is mainly due to the inclusion of transitions with low Einstein-A coefficients (<10⁻² s⁻¹). These low Einstein-A coefficient transitions are expected to have a negligible contribution to the spectrum, unless the transitions with high Einstein-A values have high optical depths of $\tau > 100$.

To identify unblended, high-energy lines or high-energy lines that are blended together, we cross-listed these lines with all the transitions available (including Einstein-A coefficients <10⁻² s⁻¹) for the molecules detected in DR Tau. As mentioned above, the transitions with low Einstein-A coefficients may have a significant contribution to the spectra if the transitions with higher Einstein-A values have high optical depths. Consequently, we also considered using these transitions when identifying the unblended lines. In this process, we only included these lines in the wavelength regions in which molecules have been detected. For C₂H₂, CO₂, and HCN this means that we cross-listed with all the available transitions in the region between 13.6 and 16.3 µm, whereas for OH we cross-listed with all transitions at wavelengths >13.6 µm. In addition, we cross-listed with all the low-energy (<6000 K) H₂O transitions. If no transitions of any of the molecules are located within five resolution elements enclosing a high-energy H₂O transition, we considered the transition to be useful for this analysis. Our cross-listing yielded a total of 24 potentially useful transitions, of which 11 are, upon visual inspection, strongly detected in both DR Tau and CI Tau. Of those 11 lines, six are blended together in three pairs and are, subsequently, treated as one.

### 3.2.2. Large disk template: CI Tau

We scaled the spectrum of CI Tau to match the flux level of the hot H₂O transitions in DR Tau by accounting for the different distances and scaling by the averaged H₂O line luminosity of the unblended lines of $L_{H_2O,DR Tau}/L_{H_2O,CI Tau} \approx 4.88$. Figure 4 displays the (scaled) spectra of DR Tau (grey) and CI Tau (black) placed at a small flux offset with respect to the residuals in red over a large portion of the 10.0–27.5 µm wavelength region. Similar to Banzatti et al. (2023a), various cold H₂O lines have excess flux present in the residuals, suggestive of a second, colder H₂O emission component, likely the result of efficient radial drift. By fitting a simple LTE slab model to the residuals, which is also displayed in Fig. 4, we find that the second component has an excitation temperature of $T = 375$ K, traces a column density of $N = 2.0 \times 10^{19}$ cm⁻², and has an emitting radius of $R_{em} = 0.91$ au.

Two lines just short of ~24 µm, indicated by the black, dashed rectangle in Fig. 4, are not well fitted by the residual slab model. These lines are also visible in the work of Banzatti et al. (2023a), who found that these are best fitted by an even colder component at $T \sim 170$ K, which most prominently emits at even larger wavelengths, ≥30 µm (Zhang et al. 2013; Blevins et al. 2016; Banzatti et al. 2023a).

## 3.3. Cold and warm H₂O reservoirs: multi-component slab models

As shown in Sect. 3.2.2, the pure rotational H₂O spectrum (≥10 µm) of DR Tau is best described by multiple, at least three, temperature components. In addition, following Sect. 3.1,

Fig. C2 shows all lines available for our analysis without a lower limit imposed on $A_{ul}$. The lines highlighted by the black circles have been used for the scaling of CI Tau.

<footnote n="4">⁴ https://zenodo.org/records/12703820: Fig. C1 shows only the lines with $A_{ul} \ge 10^{-2}$ s⁻¹ and wavelengths above 10 µm, whereas</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Spectra (across 13.4–24.0 µm) of DR Tau (grey) and CI Tau (black, scaled; see Sect. 3.2.2) shown together with the residual spectrum (in red) of DR Tau after subtraction of the scaled spectrum of CI Tau. The best fitting H₂O slab model (T = 375 K) to the residuals is shown in blue. The black dashed box just shortward of ~24.0 µm indicates the pair of lines identified by Banzatti et al. (2023a), hinting at a third component (~170 K) needed to fully explain the observed H₂O reservoir.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Cartoon visualising the different radial and vertical temperature gradients tested in this work. A top down version is shown in the top half of the figure, whereas an edge-on view is shown in the bottom half. The coloured arrows indicate the emission from the different components, with red showing the emission from the hottest component, light blue for the component with the intermediate temperature, and dark blue for the coldest component. The dashed arrows indicate emission originating from spatially shielded regions.</figcaption></figure>

the hotter (colder) components are thought to trace emission at shorter (longer) wavelengths with smaller (larger) emitting areas. A next step in the analysis of H₂O would be to combine these notions: i.e., we fitted multiple slab models

to the spectrum, consisting of components with decreasing temperatures and increasing emitting areas. The fitting is carried out using the Monte-Carlo Markov Chain implementation EMCEE (Foreman-Mackey et al. 2013), using 250 walkers and
150 000 iterations. We test three different approaches for fitting multiple components, see also Fig. 5.

In the first approach (I) we assumed that the temperature only varies radially, accounting for the decreasing temperature gradient found in planet-forming disks. We consider that each component has a separate emitting region, with hotter components emitting from the smallest regions, closest to the host star. A simple radial temperature gradient is introduced by considering a weighted sum of the components, where the weights correspond to the emitting area (A), $F_{\text{total}} = \sum_i F_i A_i$. The flux of each component, $F_i$, is determined by a temperature $T_i$ (in K) and a column density $N_i$ (in cm$^{-2}$), whereas the emitting area is parameterised by an emitting radius $R_i$ (in au). For our three components, this yields the following total flux:

$$F_{\text{total}} = F_1\pi\left(\frac{R_1}{1\text{au}}\right)^2 + F_2\pi\left[\left(\frac{R_2}{1\text{au}}\right)^2 - \left(\frac{R_1}{1\text{au}}\right)^2\right] \\ + F_3\pi\left[\left(\frac{R_3}{1\text{au}}\right)^2 - \left(\frac{R_2}{1\text{au}}\right)^2\right]. \quad (2)$$

As planet-forming disks also have a vertical temperature gradient, we allowed the colder components in the second approach (II) to also emit from the overlapping regions with the hotter components. If we are not taking shielding into account and assume that the various components have overlapping emitting areas, the total flux simply becomes:

$$F_{\text{total}} = F_1\pi\left(\frac{R_1}{1\text{au}}\right)^2 + F_2\pi\left(\frac{R_2}{1\text{au}}\right)^2 + F_3\pi\left(\frac{R_3}{1\text{au}}\right)^2. \quad (3)$$

In the final approach (III), we accounted for shielding of the colder, deeper components by the optical depth of the hotter components. In the case of shielding, the flux is attenuated by a factor of $\exp(-\sum_i \tau_i)$, where the sum is taken over the optical depths of all the components with a higher temperature. The optical depth is taken to be wavelength-dependent determined for all transitions combined using the corresponding excitation temperature and column density. As shown in Sects. 3.1 and 4.1, and in Table 1, the mutual line shielding by neighbouring H$_2$O transitions does not have a significant impact on the slab models and, hence, we only account for the spatial shielding of colder components by the hotter ones (see Fig. 5). The total flux of the this approach, in which we account for spatial line shielding of the colder components by the hotter one, becomes:

$$F_{\text{total}} = F_1\pi\left(\frac{R_1}{1\text{au}}\right)^2 + F_2\pi\left(\frac{R_1}{1\text{au}}\right)^2 \exp(-\tau_1) \\ + F_2\pi\left[\left(\frac{R_2}{1\text{au}}\right)^2 - \left(\frac{R_1}{1\text{au}}\right)^2\right] + F_3\pi\left(\frac{R_1}{1\text{au}}\right)^2 \exp(-\tau_1 + \tau_2) \\ + F_3\pi\left[\left(\frac{R_2}{1\text{au}}\right)^2 - \left(\frac{R_1}{1\text{au}}\right)^2\right] \exp(-\tau_2) \\ + F_3\pi\left[\left(\frac{R_3}{1\text{au}}\right)^2 - \left(\frac{R_2}{1\text{au}}\right)^2\right]. \quad (4)$$

For the two component version, all terms including $F_3$ can be neglected.

For all three approaches, we considered three components and imposed the following constraints: the temperature must decrease with each component, whereas the emitting radius must increase. We kept the column density as a completely free parameter. In addition, we also show a model with two components for the third method, which allows for a thorough comparison between the models. The results for the different approaches are presented in the following subsections.

<table><caption>**Table 2.** Best-fit parameters for the multiple component slab model fits, including radial and vertical temperature gradients.</caption><thead><tr><th>Component</th><th><i>T</i> (K)</th><th>log<sub>10</sub> (<i>N</i>)<sup>(<i>a</i>)</sup></th><th><i>R</i><sub>em</sub> (au)</th><th><i>N</i></th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Approach I</td></tr><tr><td>1</td><td>806<sup>+289</sup><sub>-154</sub></td><td>19.2<sup>+1.5</sup><sub>-0.5</sub></td><td>0.35±0.15</td><td>1.40 × 10<sup>45</sup></td></tr><tr><td>2</td><td>468<sup>+79</sup><sub>-85</sub></td><td>18.5<sup>+0.9</sup><sub>-0.4</sub></td><td>1.18<sup>+0.27</sup><sub>-0.18</sub></td><td>3.03 × 10<sup>45</sup></td></tr><tr><td>3</td><td>181<sup>+36</sup><sub>-43</sub></td><td>17.9<sup>+2.2</sup><sub>-0.7</sub></td><td>6.45<sup>+2.26</sup><sub>-2.64</sub></td><td>2.30 × 10<sup>46</sup></td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Approach II</td></tr><tr><td>1</td><td>796<sup>+268</sup><sub>-141</sub></td><td>19.2<sup>+1.4</sup><sub>-0.5</sub></td><td>0.36<sup>+0.14</sup><sub>-0.15</sub></td><td>1.50 × 10<sup>45</sup></td></tr><tr><td>2</td><td>466<sup>+77</sup><sub>-84</sub></td><td>18.5<sup>+0.9</sup><sub>-0.4</sub></td><td>1.14<sup>+0.24</sup><sub>-0.17</sub></td><td>2.78 × 10<sup>45</sup></td></tr><tr><td>3</td><td>184<sup>+29</sup><sub>-42</sub></td><td>17.6<sup>+2.2</sup><sub>-0.7</sub></td><td>7.98<sup>+4.49</sup><sub>-3.52</sub></td><td>1.68 × 10<sup>46</sup></td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Approach III – two components</td></tr><tr><td>1</td><td>712<sup>+130</sup><sub>-100</sub></td><td>19.1<sup>+0.8</sup><sub>-0.3</sub></td><td>0.50<sup>+0.14</sup><sub>-0.12</sub></td><td>2.02 × 10<sup>45</sup></td></tr><tr><td>2</td><td>341<sup>+74</sup><sub>-78</sub></td><td>18.5<sup>+1.0</sup><sub>-0.4</sub></td><td>1.83<sup>+0.54</sup><sub>-0.32</sub></td><td>7.97 × 10<sup>45</sup></td></tr></tbody></table>

Notes. <sup>(<i>a</i>)</sup>: <i>N</i> is given in units of cm$^{-2}$.

### 3.3.1. Radial temperature gradient

The top part of Table 2 contains median values yielded by the MCMC for the first approach. The uncertainties are given, respectively, by the 16th and 84th percentiles. For this first approach, we retrieve temperatures of $806_{-154}^{+289}$ K, $468_{-85}^{+79}$ K, and $181_{-43}^{+36}$ K for the different components. In addition, we obtain combinations for the column density (log$_{10}$ (<i>N</i>), with <i>N</i> in units of cm$^{-2}$) and emitting radius of, respectively, $19.2_{-0.5}^{+1.5}$ and $0.35_{-0.4}^{+0.15}$ au, $18.5_{-0.18}^{+0.27}$ au, and $17.9_{-0.7}^{+2.2}$ and $6.45_{-2.64}^{+2.26}$ au. The best fitting model is displayed in the top panel of Fig. 6. The corner plots are available online$^5$.

### 3.3.2. Radial and vertical temperature gradient: no shielding

The resulting values for the second approach are shown in the second part of Table 2. Compared with approach I, assuming only a radial gradient with no spatial overlap, our best fit parameters do not significantly change. In particular, the main difference can be found in the emitting radius of the third, coldest component. However, the values still agree within the given uncertainties and show a clear gradient. The final model is shown in the middle panel of Fig. 6. The corresponding corner plots are available online$^6$.

### 3.3.3. Radial and vertical temperature gradient: accounting for spatial shielding

The results for the three component and two component fits, using the third approach, are shown in respectively the third and bottom parts of Table 2. For the three component fit (also visible

<footnote>$^5$ https://zenodo.org/records/12703820: Fig. E1 displays the corner plot for approach I.</footnote>

<footnote>$^6$ https://zenodo.org/records/12703820: Fig. E2 shows the corner plot for approach II.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Multi-component slab model fits for the different methods of radial gradients: method 1 (radial temperature gradient, without overlapping regions) is shown in the top panel, whereas the methods 2 (radial gradients with overlapping regions, but no line shielding) and 3 (radial and vertical temperature gradients) are displayed in the middle and bottom panels, respectively. The full model is in each panel shown in green, whereas the components are coloured depending on the corresponding temperatures. Red indicates the hottest component, followed by the light blue and dark blue ones. The green horizontal bars indicate the regions used in the $\chi^2_{red}$-fits.</figcaption></figure>

in the bottom panel of Fig. 6), we find that the best fit parameters are similar to those of the second approach. That means that the main difference can once more be found in the emitting

radius of the third component. It must be noted that the values for all three approaches agree with one another within the given uncertainties. In accordance with Pontoppidan et al. (2024), we
<figure><img src="image_7.png" /></figure>

**Fig. 7.** Cartoon visualising the expected emission regions, based on the derived excitation temperatures, emitting radii, and optical depths, of the multiple components of H₂O. Shown are also the expected emission locations of CO, CO₂, HCN, and C₂H₂ adopted from Temmink et al. (2024). The CO is decomposed in two components: one broad component (BC) tracing the Keplerian rotation of the disk, and a narrow component (NC) tracing a disk wind.

find that the results for this approach can be well described by the power-law $T(R_{em}) \sim 500(\frac{R_{em}}{1\ \text{au}})^{-0.5}$ K for the emitting layer. For the two component fit, we see that the temperature of the first (hottest) component and the second component are lower by ≥100 K. It is clear that the coldest component, $T \sim 180$ K, is not considered by the two component fit. A figure displaying the two component fit is available online<fnref n="7" />. The column density and the emitting radius of the first component have, respectively, slightly decreased and increased with respect to the three component fit. On the other hand, only the emitting radius of the second component increased from ~1.20 au to ~1.83 au. The corner plots for the three and two components fits are available online<fnref n="8" />.

We find that all models underproduce the flux for a range of H₂O transitions between ~14 µm and ~18 µm. This region is dominated by both the first (hottest) and second (warm) component, suggesting that either one or both components do not properly describe this region and slightly different values for the parameters would yield a better fit. Nonetheless, the fit provides an overall good description of the observed rotational spectrum of H₂O in DR Tau.

Finally, Fig. 7 displays an cartoon of the disk, indicating the expected emission locations from our multi-component fit. We have also included the expected locations from CO, CO₂, HCN, C₂H₂ as found in Temmink et al. (2024) (see also their Fig. 12), based on their excitation temperatures and emitting areas.

### 3.4. Line pair ratios

The inferred column density of the slab models can be tested and further constrained by comparing the peak fluxes of H₂O line pairs which have the same value for $E_{up}$, but different values for $A_{ul}$ (Gasman et al. 2023). As the strength of these lines and, therefore, their ratio depends primarily on the line opacity, they can be used to approximate as independent constraints on the column density and provide a sanity check for the column densities retrieved with the slab models. If both lines are optically thin, the flux ratio is expected to converge to the ratio of the respective values for $A_{ul}$, whereas the flux ratio will deviate from

<figure><img src="image_8.png" /></figure>

**Fig. 8.** Logarithm of the line flux ratios as a function of the logarithm of the column densities (in cm⁻²) for the line pairs with similar upper level energies. The ratios for the slab models with $T=300$ K are shown in blue, whereas those for $T=500$ K and $T=900$ are shown in green and red, respectively. The horizontal black line indicates the line ratio taken from the continuum-subtracted spectrum of DR Tau. The dashed, vertical lines indicate where the different model ratios equal the observed line ratio. In addition, the dotted lines and the column densities listed in parenthesis indicate the results for the top line pair after subtracting off the contribution of the component with the lowest temperature (see Sect. 3.3).

the $A_{ul}$ ratio if one of the lines becomes optically thick. We use the same, pure rotational line pairs as Gasman et al. (2023) (see their Table 2) to investigate the H₂O column density. Some of the lines are part of line clusters and are, subsequently, blended, which means that it cannot be concluded from the flux ratio alone which of the lines are optically thick, if the flux ratio deviates from that of $A_{ul}$.

Figure 8 shows the model line ratios for different temperatures (300, 500, and 900 K) and column densities ($14 \le \log_{10}(N) \le 22$, with N in units of cm⁻²). The ratios become constant for column densities of $\le 10^{17}$ cm⁻². For these column densities we can expect both lines to be optically thin and, subsequently, their ratios to be equal to the values for $A_{ul}$. At larger column densities (at least) one of the lines becomes optically thick and the ratio deviates from that for $A_{ul}$. For even larger column densities ($\ge 10^{20}$ cm⁻²) the ratio becomes constant again, the effect of both lines being fully saturated.

The black horizontal line shown in Fig. 8 is the line ratio obtained from the JWST-MIRI spectrum. The coloured, dotted vertical lines are the intercepts between the black horizontal line and the model ratio curves. Based on the intercept, we can expect the actual H₂O column density to take on values of $\le 10^{19.4}$ cm⁻², $\le 10^{18.4}$ cm⁻², and/or $\le 10^{18.0}$ cm⁻² for the different temperatures ($T=300, 500$, and $900$ K), respectively. Using the results from the multi-component fits, we can infer the contributions from

<footnote n="7">⁷ https://zenodo.org/records/12703820: Fig. D1 shows a comparison between the three and two component fits. The lack of the coldest component is visible in this figure, as the lines just short of 24 µm are not well fitted by the model.</footnote>

<footnote n="8">⁸ https://zenodo.org/records/12703820: Figs. E3 and E4 display the corner plots for the three and two component fit of approach III, respectively.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Zoom-ins on the isolated H₂ ¹⁸O lines, highlighted by the red vertical line. A simple H₂ ¹⁸O slab model (T = 350 K, log₁₀ (N) = 17.66, and R<sub>em</sub> = 1.0 au) is shown in red, whereas the best models of H₂O (three component fit accounting for both the radial and vertical temperature gradient) and OH are shown in blue and magenta, respectively. For each transition, we have listed the upper level energies and Einstein-A coefficients.</figcaption></figure>

each component to the lines used for both line pairs. For the bottom pair (11<sub>5,6</sub>–10<sub>4,7</sub>/11<sub>5,6</sub>–10<sub>2,9</sub>), both transitions consist of contributions from the hotter two components, so the *T* ~ 300 K models can be ignored. For both transitions, the lower temperature (*T* ~ 475 K) has a slightly larger contribution to the flux. Based on this line pair and the temperature from both contributing components (*T* ~ 830 K and *T* ~ 475 K), we infer that the column density of the combination of these two H₂O components must be log₁₀ (*N*) ≤ 18.5 within ~1.2 au. For the top pair (8<sub>7,2</sub>–7<sub>4,3</sub>/8<sub>7,2</sub>–7<sub>6,1</sub>), however, the lowest temperature component (*T* ~ 185 K) also has a significant contribution to the line flux of one of the lines. As only one of these lines has a contribution from the lowest temperature component, this pair does not provide strong constraints. After subtracting off the model corresponding to this temperature, we do infer a slightly more stringent constraint on the H₂O column density (in cm⁻²) of log₁₀ (*N*) ≤ 18.4 for the first two components combined, compared to value derived for the other line pair. Even though H₂O is optically thick, these line pair ratios provide good constraints on the column densities for the first two components. We note that the derived limit agrees well with those obtained for the second, warm component, whereas it falls just outside the uncertainties for the hottest, most optically thick component.

### 3.5. Searching for H₂ ¹⁸O

We are inconclusive about the detection of H₂ ¹⁸O in DR Tau. Some of the transitions, in particular those at 25.65, 26.83, and 26.99 µm, look promising (see Fig. 9). Accounting for the large error in subband 4C (see Table A.1), the peak fluxes suggest maximum detection levels of ≲4σ and should be statistically significant for the confirmation of a detection. However, the lack of flux seen for the transitions at ~22.03 µm and ~23.17 µm provide

the strongest argument against a detection of H₂ ¹⁸O in DR Tau. Re-observing DR Tau with JWST-MIRI with the aim of achieving higher signal-to-noise ratios across Channel 4 may improve our chances of detecting H₂ ¹⁸O.

Using our three component models, we find that the majority of the H₂O emission features are optically thick. While transitions with low Einstein-A values may be optically thin and can be used to derive better constraints on the total H₂O vapour reservoir present in the inner disk, we note that these lines are very scarce in the spectrum of DR Tau. In total, we identify four clean emission features (between 12.496–12.506, 14.377–14.392, 16.362–16.375, and 19.370–19.392 µm) that are optically thin (τ <1) according to our models and not blended with other (unidentified) emission features. All these transitions trace mainly the hottest component. Subsequently, a fit to these transitions only may provide a better constraint on the column density for the hottest component and observations of the optically thin isotopologues, such as H₂ ¹⁸O, could help to infer a better estimate of the column density and the total number of H₂O molecules present in the inner disk.

Calahan et al. (2022) list eight lines in the 19.0–27.0 µm wavelength range that may be used for the detection of H₂ ¹⁸O in JWST-MIRI spectra (see their Table 1). We, however, note that the line at 20.01 µm is blended with OH and can likely not be used for the detection of H₂ ¹⁸O in many T-Tauri disks. On the other hand, through slab model exploration we discovered a blended pair of H₂ ¹⁸O lines that can be used instead. The line pair (8<sub>8,0</sub>–7<sub>7,1</sub> and 8<sub>8,1</sub>–7<sub>7,0</sub>) is located at 25.65 µm and have the same upper level energies (*E*<sub>up</sub> ~2544.9 K) and Einstein-A coefficients (*A*<sub>ul</sub> ~29.08 s⁻¹), suggesting that both lines contribute half of the observed line flux. The individual lines are displayed in Fig. 9. The H₂ ¹⁸O emission is highlighted by a simple H₂ ¹⁸O slab model (T = 350 K, log₁₀ (*N*) = 17.36, and R<sub>em</sub> = 1.0 au),
where we have accounted for the isotopologue ratio in the local interstellar medium (ISM) of ¹⁶O/¹⁸O~550 (Wilson 1999) in the column density, i.e. we have divided the found column density of H₂O in the wavelength region of 23.0–27.0 µm (log (N) ~ 20.1, see Table 1) by 550. We show our best fitting H₂ ¹⁶O and OH slab models (see Table 1 and Sect. 3.6.1) for further identification. A figure displaying the full 17.0–27.5 µm wavelength region, including the aforementioned H₂ ¹⁶O, H₂ ¹⁸O, and OH slab models, is available online<fnref n="9" />.

Even though we are inconclusive about the detection of H₂ ¹⁸O, the used slab model yields underproduced fluxes for the most promising peaks. This may suggest that, for the used temperature, the input column density is not high enough. As the column density was set to account for the isotopologue ratio (¹⁶O/¹⁸O~550), this indicates optically thick H₂O emission, which is in agreement with our models. We note that an underestimate of the column density for H₂ ¹⁸O may also point to a deviation of the isotopologue from the ISM value. As discussed in Calahan et al. (2022), the gas may be enriched in H₂ ¹⁸O in the disk layers where CO becomes optically thick and self-shielding against photodissociation becomes important. As the main isotopologue, C¹⁶O, is more abundant than C¹⁸O, there will be a layer where C¹⁶O will be self-shielded, whereas C¹⁸O is still being dissociated. The gas in this layer will in that case be enriched in ¹⁸O atoms, which may end up in H₂O-molecules, enhancing the H₂ ¹⁸O column. Additionally, UV-shielding by the main isotopologue, H₂ ¹⁶O, can also enhance the H₂ ¹⁸O abundance.

## 3.6. Other emission features

Besides emission from H₂O, the spectrum of DR Tau is rather molecule rich. CO, CO₂, HCN, and C₂H₂ have been thoroughly discussed been in Temmink et al. (2024), so the upcoming subsections discuss the other molecular features potentially present in the spectrum, including OH and the larger hydrocarbons that have been observed in the disks around very low-mass stars (e.g. Tabone et al. 2023).

### 3.6.1. OH

Since we do not detect any OH lines at wavelengths <13 µm, similar to Gasman et al. (2023), we do not expect prompt emission, following H₂O photodissociation (Najita et al. 2010), to play a role in the observed OH emission in DR Tau. We fit the OH emission, using a reduced χ²-minimisation (see also Sect. 3.1), over the entire wavelength range of 10 µm to 27.5 µm. We have used the same grid for T, log₁₀(N), and Rₐm as in Sect. 3, except we extended the grid for the excitation temperature up to 4000 K. The emission was fit to the residual spectrum after subtracting the resulting H₂O model spectrum taking both a radial and vertical temperature gradient into account (see Sect. 3.3.3). For OH, we retrieve a temperature of T=1875[775] K, a column density (in cm⁻²) of log (N)=13.5[3.3], and emitting radius of Rₐm=9.49[0.52] au. As the emission is optically thin (τ ≤ 0.005) for this combination of excitation temperature and column density, the emitting radius is not well constrained, as indicated by the large lower uncertainty. We note that the upper uncertainty of the emitting radius and, therefore, the lower constraints on the column density are set by the limits of our

grids. For these parameters, we derive that the number of OH molecules present in the inner disk of DR Tau is N ~ 2.0 × 10⁴². Figures displaying parts of the best fitting model, including the wavelength regions used in the fitting, are available online<fnref n="10" />. Similar to Gasman et al. (2023), the OH lines probe larger gas temperatures compared to H₂O, but lower column densities. As the excitation of the OH lines includes prompt emission (Tabone et al. 2021, 2024) and collisional excitation or chemical pumping (Zannese et al. 2024) through the O+H₂→OH+H reaction, it is unlikely that the observed transitions follow a single excitation temperature. Therefore, the inferred temperature from the slab model is not related to a kinetic temperature. While they cannot infer a kinetic temperature, we note that slab models are a powerful tool for the identification of OH transitions in the JWST-MIRI/MRS spectra of disks. Detailed models, including the different excitation pathways, can be used to further explore the excitation properties and to infer other key disk parameters (see, for example, Tabone et al. 2024).

### 3.6.2. Atomic and molecular hydrogen

We list the detection of three atomic hydrogen (H) transitions (10–6, 6–5, and 8–6) and the potential detection of four molecular hydrogen (H₂) transitions: 0, 0 S(1), S(2), S(3), and S(4). As the H₂ lines are either blended with H₂O or potentially blended with unidentified features, their detections are not fully secure. Higher transitions (e.g. S(5), S(6), S(7), and S(8)) may also be present, however, their detections are even less certain due to stronger blending with H₂O and CO emission features. Figures displaying the (potentially) observed transitions of H and H₂ (S(1), S(2), S(3), and S(4)) are available online<fnref n="11" />. We also see no evidence of extended H₂ emission coming from the DR Tau disk in the form of a disk wind or an outflow, however, we do observe extended H₂ emission (most prominent in the S(1) transition) from the cloud surrounding DR Tau, as was previously observed by Thi et al. (2001) in millimetre single dish CO J = 3–2 transitions offset by ~2 km s⁻¹. A figure visualising the background emission is available online<fnref n="12" />. Visual inspection of the spectral cube shows that the flux of the cloud emission (FₐS(1) ~ 2 × 10² MJy sr⁻¹) is a factor ~100 lower than the flux at the position of DR Tau (FₐS(1) ~ 2 × 10⁴ MJy sr⁻¹). As we approximate our background using an annulus, the H₂ emission from the cloud will be captured in that annulus. However, as the flux is significantly lower than the source flux, the background subtraction does not impact our detections.

Using the residual flux, after subtracting the corresponding H₂O slab model fit, we determined an upper limit on the column density of H₂ based on the 0, 0 S(4) transition (Eₐup. = 3474.5 K). Assuming that H₂ comes from the inner most region, we adopt a temperature of T = 1000 K and an emitting radius of Rₐem = 0.20 au, as this also allows us to make direct comparisons with the total number of molecules of both CO and H₂O (as obtained in regions 1 and 2 for the ro-vibrational transitions, see Sect. 3.1) derived in these innermost regions. To calculate the column density under the assumption of optically thin emission, we follow

<footnote n="9">⁹ https://zenodo.org/records/12703820: Fig. F1 displays the H₂¹⁶O, H₂¹⁸O, and OH slab models over the 17.0–27.5 µm wavelength region.</footnote>

<footnote n="10">¹⁰ https://zenodo.org/records/12703820: Figs. B1-4 to B1-6 display parts of the best fitting OH slab model.</footnote>

<footnote n="11">¹¹ https://zenodo.org/records/12703820: Figs. G1 and G2 display, respectively, the (potentially) observed transitions of H and H₂.</footnote>

<footnote n="12">¹² https://zenodo.org/records/12703820: Fig. G3 displays the H₂ emission coming from the background cloud.</footnote>
the description given in Goldsmith & Langer (1999):

$$ \frac{4\pi F}{hcv\Omega g_{up}A_{ul}} = \frac{N_{tot}}{Q(T)} \exp\left(-\frac{E_{up}}{k_{B}T}\right), \qquad (5) $$

where *F* is the integrated flux (in erg s⁻¹ cm⁻²), *ν* the frequency (in cm⁻¹), *g*<sub>up</sub> the upper level degeneracy, *Q*(*T*) the partition function at a temperature *T*, and *E*<sub>up</sub> the upper level energy (in K). Using the resolution element at the location of the S(4) transition and the two resolution elements directly next to it, we derive *F*=3.9×10⁻¹⁵ erg s⁻¹ cm⁻². As mentioned above, the H₂ flux might be slightly affected by the annulus background subtraction due to cloud emission. However, as the background emission is significantly lower than the on-source emission, we expect the derived flux to be minimally impacted. This yields a column density of *N*<sub>tot</sub> ~ 7.6 × 10²⁵ cm⁻² and a total number of molecules, assuming an emitting radius of 0.20 au, of *N*<sub>H₂</sub> ~ 2.1×10⁵¹. We note once more that these derived column densities and total number of molecules hold for the upper layers of the disk, above the region where the dust emission becomes optically thick, or, in the case of optically thick lines, above the height where the emission lines become optically thick (see, for example, Bosman et al. 2022).

### 3.6.3. Non-detections

In an attempt to identify as many molecular features as possible, we have also looked for emission signatures of various other molecules and isotopologues. These species are comprised of ¹³CO₂, CH₄, NH₃, CS, H₂S, SO₂, and various larger hydrocarbons, including ¹³CCH₂. We have obtained the required spectroscopic data for these species through the HITRAN database (Gordon et al. 2022). Table 3 summarises the parameters (assuming an emitting radius of *R*<sub>em</sub> = 1.0 au) of simple LTE slab models with fixed temperature and column density that have been used to search for these species. In addition, the table lists the best fit parameters for CO₂, HCN, and C₂H₂, as found in Temmink et al. (2024). Figures displaying the slab models are available online<fnref n="13" />.

We are not able to confidently detect any of the listed molecular species. The non-detection of some the molecules, for example, ¹³CO₂, CH₄, SO₂, and the larger hydrocarbons, becomes apparent from the slab models, mostly due to the lack of a molecular continuum. For the other molecules, NH₃, CS, H₂S, and ¹³CCH₂, a (non-)detection cannot be confirmed, largely due to leftover broad features (especially visible at the location of the silicate feature at ~10 µm), which are potentially of molecular nature, in the spectrum. We also do not detect the larger hydrocarbons, mostly due to the lack of molecular continua. While the slab model<fnref n="14" /> of C₆H₆ looks promising, the presence of C₆H₆ cannot be confirmed in the spectrum of DR Tau, as the lack of emission observed for the other hydrocarbons makes the presence of C₆H₆ less likely.

For most of the non-detected molecules, the upper limit on the total number of molecules is on the order of *N* ~ 10⁴²–10⁴³. Only the upper limit for H₂S is larger (*N* ~ 7.0×10⁴⁴), whereas that for C₆H₆ is significantly lower (*N* ~ 1.3×10³⁸). Using the number of molecules from our three component H₂O

<table><caption>**Table 3.** Slab models parameters of the various detected (CO₂, HCN, C₂H₂; Temmink et al. 2024) and non-detected species, using *R*<sub>em</sub>=1 au.</caption><thead><tr><th>Mol.</th><th>Region (µm)</th><th>T (K)</th><th>log<sub>10</sub> (N)<sup>(α)</sup></th><th>N</th></tr></thead><tbody><tr><td>¹²CO₂<sup>(β)</sup></td><td>13.60–16.30</td><td>325<sup>+50</sup><sub>-100</sub></td><td>17.4<sup>+0.7</sup><sub>-0.2</sub></td><td>4.96 × 10<sup>43</sup></td></tr><tr><td>HCN<sup>β</sup></td><td>13.60–16.30</td><td>900<sup>+50</sup><sub>-50</sub></td><td>14.7<sup>+2.1</sup><sub>-0.6</sub></td><td>8.22 × 10<sup>42</sup></td></tr><tr><td>C₂H₂<sup>(β)</sup></td><td>13.60–16.30</td><td>775<sup>+150</sup><sub>-625</sub></td><td>15.0<sup>+5.6</sup><sub>-1.7</sub></td><td>1.48 × 10<sup>42</sup></td></tr><tr><td>¹³CO₂</td><td>15.30–15.50</td><td>300</td><td>15.57</td><td>≤2.6 × 10<sup>42</sup></td></tr><tr><td>CH₄</td><td>7.60–7.75</td><td>700</td><td>15.75</td><td>≤4.0 × 10<sup>42</sup></td></tr><tr><td>NH₃</td><td>8.50–9.50</td><td>500</td><td>15.50</td><td>≤2.2 × 10<sup>42</sup></td></tr><tr><td>CS</td><td>7.50–8.30</td><td>500</td><td>15.50</td><td>≤2.2 × 10<sup>42</sup></td></tr><tr><td>H₂S</td><td>7.00–8.00</td><td>500</td><td>18.00</td><td>≤7.0 × 10<sup>44</sup></td></tr><tr><td>SO₂</td><td>7.10–7.62</td><td>500</td><td>15.50</td><td>≤2.2 × 10<sup>42</sup></td></tr><tr><td>¹³CCH₂</td><td>13.60–13.85</td><td>750</td><td>15.00</td><td>≤7.0 × 10<sup>41</sup></td></tr><tr><td>C₂H₄</td><td>10.50–10.60</td><td>700</td><td>16.00</td><td>≤7.0 × 10<sup>42</sup></td></tr><tr><td>C₂H₆</td><td>12.00–12.34</td><td>700</td><td>17.00</td><td>≤7.0 × 10<sup>43</sup></td></tr><tr><td>C₄H₂</td><td>15.80–16.00</td><td>700</td><td>15.50</td><td>≤2.2 × 10<sup>42</sup></td></tr><tr><td>C₆H₆</td><td>14.70–15.00</td><td>700</td><td>11.25</td><td>≤1.3 × 10<sup>38</sup></td></tr></tbody></table>

**Notes.** (α): *N* is given in units of cm⁻². (β): The listed values for CO₂, HCN, and C₂H₂ are taken from Temmink et al. (2024) and are shown here for completion.

fits (*N* ~ 10⁴⁵–10⁴⁶), most of the non-detected molecules are a factor ≥10²–10⁴ less abundant than H₂O. H₂S and C₆H₆ are factors of, respectively, ≥10–100 and ≥10⁷–10⁸ less abundant. On the other hand, the upper limits are of a similar level compared to the total number of molecules derived for CO₂, HCN, and C₂H₂. Once again, the upper limit for H₂S is slightly larger, whereas that of C₆H₆ is a few orders of magnitude lower. For comparison, the abundance of CO₂ is found to be a factor ~60 lower than the second, warm component of H₂O.

## 4. Discussion

### 4.1. The need for H₂O line overlap

For slab model fits of molecules with a well defined Q-branch (see Tabone et al. 2023), line overlap (i.e. mutual shielding of lines) is often included to properly fit the spectra. The spectrum of H₂O, on the other hand, consists of well defined, separate/isolated transitions. In the disk of AS 209, Muñoz-Romero et al. (2024) also tested the need for line overlap of multiple molecular species. They conclude that a proper treatment of the line overlap is necessary when analysing the observed molecular species together.

In Sect. 3.1, we test the need for mutual shielding by the H₂O transitions alone. We test this by creating slab models without and with mutual line shielding. As can be seen in Table 1, the differences between the slab models without (top part of the table) and those with (bottom part) mutual line overlap are negligible when fitting H₂O alone. The largest differences are found for temperatures derived for the regions spanning the ro-vibrational transitions, which is likely due to the lines being more densely packed in this part of the spectrum compared to those in the pure rotational part. We note that the residuals for the models without and with mutual line shielding are quite similar. We find that one model will provide a better fit for some transitions, whereas other transitions are better fitted by the other model. Therefore, we

<footnote n="13">¹³ https://zenodo.org/records/12703820: Figs. G4-1 to G4-6 and Figs. G5-1 to G5-5 display, respectively, the slab models of the listed species and the larger hydrocarbons.</footnote>

<footnote n="14">¹⁴ https://zenodo.org/records/12703820: Fig. G5-5 displays the C₆H₆ slab model.</footnote>
conclude that the inclusion of mutual shielding of neighbouring H₂O lines is not of particular importance for the pure rotational transitions of H₂O, when fitting them alone, and that the slab models without mutual line shielding can be used to infer constraints on the excitation conditions. As our slab models yield different results for the rovibrational part of the spectrum, the inclusion of mutual shielding may be more important for these more densely packed transitions. Mutual line shielding is of particular importance for molecules with a well defined *Q*-branch, such as CO₂, HCN, and C₂H₂.

## 4.2. The need for a radial temperature gradient

As expected, the various methods all indicate a radial temperature gradient visible in the spectrum. Fitting, however, the H₂O spectrum across multiple, independent wavelength regions does not yield any information on the low temperature (T ~ 170–200 K, see Table 1) component, even though the fitted regions include the lines best suitable for finding this component. Subsequently, to properly investigate the various components, one could use the large disk template method implemented by Banzatti et al. (2023a) or the multi-component slab model method introduced here.

Rather surprisingly, the current approaches of the multi-component slab model fits yield the same results for all three methods: temperatures of ~800 K, ~470 K, and ~180 K, column densities of ~10¹⁹.² cm⁻², ~10¹⁸.⁵ cm⁻², and ~10¹⁷.⁴ cm⁻², and emitting radii of ~0.3 au, ~1.2 au, and ~6.5–8.1 au. For comparison, the midplane H₂O snowline in DR Tau can be estimated from a simple power law (Chiang & Goldreich 1997; Dullemond et al. 2001; van der Marel et al. 2021) involving the total luminosity ($L_{\text{tot}}=L_{\text{s}}+L_{\text{acc}} \approx 0.63 + 0.58 L_{\odot}$; Long et al. 2019; Banzatti et al. 2020), providing a radius $R_{\text{H}_2\text{O}} \sim 0.75$ au. The calculations of Mulders et al. (2015) (see their Fig. 1) suggest that the snowline is located at a radial distance of ~1–2 au. Radiative transfer modelling of the dust disks may provide further constraints on the location of the H₂O snowline. As the third component at ~180 K is at larger distances than the midplane snow surface, this suggests a curved rather than vertical snow surface. We find that the decreasing temperature profile can be well described by the powerlaw $T(\text{R}_{\text{em}}) \sim 500(\text{R}_{\text{em}}/\text{au})^{-0.5}$ K. Only the emitting radius for the final component differs between approach 1 and approaches 2 and 3, but still agrees well within the given uncertainties. The resulting temperatures agree with those given by the large scale template method, in particular the third component (T ~ 180 K) corresponds well with the one (T ~ 170 K) that Banzatti et al. (2023a) needed to explain the lines just shortward of 24 µm. To extend upon this notion, a figure showing a comparison between the multi-component fits including three and two components is available online<fnref n="15" />. The main difference is most clearly visible in the region of 23.50–24.25 µm, which includes the aforementioned cold lines. The two component model is not able to properly fit this line cluster, whereas the three component model can.

Compared to results from fitting H₂O across the different wavelength regions (see Sect. 3.1), we note that the derived column densities for these methods do not show the same behaviour. For the multi-component fit, the column densities decrease with lower temperatures, whereas they increase with temperature for

the different regions. The difference in behaviour may have arisen from combining the contribution of every component in the multi-component fits, whereas contributions of other components and/or wavelength regions have not been accounted for in the fits to the different wavelength regions.

A multi-component fit could, instead of a large disk model, also be used to investigate the enhanced reservoir of cold H₂O through radial drift. We expect that a two component fit is sufficient to fit the spectra without enhancement of cold H₂O through radial drift, whereas three component fits are required to properly fit those with the enhancement. Further work examining a large sample of disks, small and large, structured and smooth, is necessary to fully explore this thought.

The reason why all three models yield the same results is due to line optical depths. These optical depths are too high (reaching values up to $\tau \sim 3000$ for the innermost region and $\tau \sim 350$ for the first annulus for the brightest lines) for any of the shielded regions to make a significant contribution to the total model. In other words, considering only the innermost region, the optical depth of the hottest components is so high that the contributions (i.e. those with an exponential in Eq. (4)) of the lower temperature components are negligible. As discussed in Sect. 3.4, the column density (in cm⁻²) of the warmer two components (components 1 and 2) can be expected to be on the order of log(N) ~ 18.4.

When using simple LTE slab models to describe the full rotational spectrum of H₂O, it is clear that a vertical temperature gradient does not need to be included. Due to the negligible contributions from the overlapping region, a simple radial gradient (using two or three components) alone should be sufficient to represent the emitting layer in the upper region of the disk. For more sophisticated thermochemical modelling methods, more realistic column densities may be probed and the inclusion of vertical temperature gradients may become of significant importance.

## 4.3. The importance of radial drift

We find that our introduced multi-component models, with various levels of complexity, are able to yield similar results as those obtained with the method of Banzatti et al. (2023a). We note that at least three components are required to fully describe the pure rotational H₂O of DR Tau (see also Sect. 4.2). In particular, these three components, including a cold component of T ~ 180 K, are necessary to describe the quartet of H₂O lines just shortward of ~24 µm. This cold component has been suggested by Banzatti et al. (2023a) to trace an additional H₂O reservoir near the snowline, following the inward drift of icy pebbles, the subsequent sublimation, and the diffusion of H₂O vapour. Modelling works (e.g. Bosman et al. 2018; Kalyaan et al. 2021) have shown that the inward drift of solids can increase the H₂O abundance (or vapour mass) inside the snowline up to an order of magnitude during the first few million years. As our models require this cold component and the disk around DR Tau is compact in the millimetre dust (≤60 au), our multi-component analysis can be used to further study the importance of radial drift in setting the potentially enhanced observable H₂O reservoirs in the inner regions of planet-forming disks. An upcoming paper will utilise this method on a larger sample of compact disks.

## 4.4. CO versus H₂O

With both CO (see Temmink et al. 2024) and ro-vibrational H₂O (Regions 1 and 2 in Table 1) analysed, we can compare

<footnote n="15">¹⁵ https://zenodo.org/records/12703820: Fig. D1 shows a comparison between the three and two component fits. The lack of the coldest component is visible in this figure, as the lines just short of 24 µm are not well fitted by the model.</footnote>
their respective total number of molecules, assuming the two molecules are co-spatial down to the layer where the ~6 µm dust continuum becomes optically thick. The assumption that the two molecules may be co-spatial is supported by the high spectral resolution study of Banzatti et al. (2023b), where the similar line profiles suggest that the ro-vibrational H₂O lines may originate from the same region as the broad component of CO, which is thought to trace the Keplerian rotation of the disk. Using the optically thin C¹⁸O emission seen in the VLT-CRIRES observations of DR Tau, Temmink et al. (2024) derived a total number of molecules of NCO=4.1×10⁴⁴. For the C¹⁸O, they retrieved an excitation temperature of T ~975 K and an emitting radius of Rem = 0.23 au, which agree well with the best fit parameters found for the ro-vibrational H₂O transitions at the shortest wavelengths (5.0–6.5 µm, see Table 1), suggesting that the probed emission may indeed come from the same region in the disk. For this wavelength range (5.0–6.5 µm), we retrieve a total number of molecules of NH₂O ~ 7×10⁴³. These values yield a ratio of NH₂O/NCO ~ 0.17.

Taking typical ISM abundances of [O]=5.8×10⁻⁴ and
[C]=3.0×10⁻⁴, and taking into account the amounts of oxygen
and carbon that are locked up in solid form, the maximum
H₂O/CO ratio is ~1.4–2.0, if all volatile oxygen not locked up
in CO is driven into H₂O (van Dishoeck et al. 2021). This
ratio is significantly higher than we retrieve for DR Tau, sug-
gesting that the total number of H₂O molecules may be low,
instead of enhanced as expected for compact disks, in the inner
disk. However, as discussed at previous points throughout this
work (see Sects. 3.5 and 4.2), many of the observed H₂O tran-
sitions are optically thick. Subsequently, our derived ratio of
N*H₂O*/N*CO* ~ 0.16 acts as a lower limit. A factor ~10 is required
to meet the aforementioned expected ratio of H₂O/CO~1.4–2.0.
Perhaps the H₂O emission becomes optically thick well above
the layer where the dust becomes optically thick, suggesting that
we may probe less deep into disk compared to both CO (probed
with the optically thin C¹⁸O emission) and H₂. Additionally,
H₂O may self-shield against the photodissociating UV-photons
before C¹⁸O (Calahan et al. 2022), leaving a thin emitting layer
where both molecules are co-located. Higher signal-to-noise
detections of one of the rarer H₂O isotopologue (e.g. H₂ ¹⁸O)
or detections of the optically thin, low Einstein-A H₂ ¹⁶O transi-
tions (see Sect. 3.5), which probe deeper into the disk (Bosman
et al. 2022), are needed to further explore the abundance ratio of
CO and H₂O.

We note that the total number of molecules for both H₂O and
CO appear to be low with respect to that of molecular hydrogen,
H₂. For CO the ISM abundance is of the order of CO/H₂~10⁻⁴
(Bergin & Williams 2017; Lacy et al. 2017), whereas that of H₂O
is H₂O/H₂~4×10⁻⁴ (van Dishoeck et al. 2021), if most volatile
oxygen is used in the formation of H₂O. Assuming a tempera-
ture of T ~ 1000 K and an emitting radius of R<sub>em</sub> = 0.20 au,
we derive an upper limit on the total number of molecules of
N<sub>H₂</sub> ≤ 2.1×10<sup>51</sup>. This number is a factor of, respectively, 5.1×10<sup>6</sup>
and 3.0×10<sup>7</sup> higher than the total number of molecules derived
for CO and H₂O, under the assumption that all molecules probe
the same region of the disk. However, caution is warranted as H₂
is not strongly detected and we emphasise that derived number
of molecules must be treated as an lower limit.

**5. Conclusions and summary**

In this work, we thoroughly analysed the H₂O excitation tem-
peratures and column densities in the inner region of the
planet-forming disk of DR Tau, as seen with JWST-MIRI/MRS.

We have used techniques presented in previous works, as well
as tried, for the first time, a multi-component approach to fully
describe the rotational H₂O spectrum (≥10 µm). In addition, to
the H₂O analysis, we have also analysed the OH, atomic and
molecular hydrogen emission, and looked for other molecular
species. We summarise our conclusions as follows:

- By analysing the H₂O over various wavelength regions, we find that the excitation temperature decreases with wavelength, whereas the emitting radius increases. This is a clear sign of a radial temperature gradient visible in the spectrum. The need for a temperature gradient, or multiple components, is further confirmed by using the spectrum of a large disk (CI Tau) as comparison template for the spectrum of our compact disk. Various cold emission lines are not reproduced by this large disk template and are explained by slab models with low temperatures of T ~ 400 K and T ~ 200 K;

- To further test the need for multiple components, we used an MCMC approach to fit two and three components at the same time. Our methods include a simple radial gradient, a radial gradient with spatial overlap, and a radial and vertical gradient. The different methods yield no significant differences in the excitation properties (temperatures, column densities, and emitting radii) of the different components, likely due to the high optical depth of the hotter components in the overlapping regions. Consequently, we can state that simple models (i.e. models only considering a radial temperature gradient) are able to provide a proper description of the rotational H₂O spectrum and are supported by the more complex models (i.e. those including a vertical gradient). Only for disks where the hotter components are not optically thick, the results will differ and a more complex model must be used;

- We find that the temperature for the third approach can be described roughly as $T(R_{\text{em}}) \sim 500 \left(\frac{R_{\text{em}}}{1 \text{ au}}\right)^{-0.5}$ K. For DR Tau, we emphasise the need for a minimum of three components ($T_1 \sim 800$ K, $T_2 \sim 470$ K, and $T_3 \sim 180$ K), as the two component model is unable to fit the transitions at the longest wavelengths, especially those around $\sim 24$ µm. While the need for the third component around the H₂O snowline is in agreement with previous works, it is not yet possible to conclude whether the H₂O abundance is enhanced due to radial drift. A future study focusing on a consistent analysis of a larger sample of compact disks will provide more insights in the importance of radial drift;

- While OH and various atomic and molecular hydrogen transitions are detected, we are inconclusive about the detection of the isotopologue H₂ ¹⁸O, mostly due to the (current) large uncertainties at the longest wavelengths. We also do not report the detections of various other species, including: CH₄, NH₃, CS, H₂, SO₂, and larger hydrocarbons. The derived upper limits on the column densities of these non-detected species are almost all significantly lower, by factors ~10²–10⁴, than those found for H₂O in the different wavelength regions (see Sect. 3.1);

- Finally, we have compared the derived number of H₂O molecules with that of CO, derived from the optically thin C¹⁸O emission. The comparison yields a low ratio of ~0.17, which, for the hot H₂O component T ~1000 K, may point to H₂O having a relative low total number of molecules in the layer that we are able to trace. However, the strongest H₂O lines are found to be optically thick, suggesting that we may be probing less deep into the disk compared to C¹⁸O, and so this ratio must be treated as a lower limit. A detection of one
of the rarer isotopologues (e.g. H₂ ¹⁸O) is required to better constrain their relative number of molecules down to the layer where the dust continuum becomes optically thick.

This work has shown the need of fitting multiple components to explain the rotational spectrum of H₂O. Future works focusing on the multiple-component analysis of H₂O in larger samples of disks can help us further understand the H₂O reservoir present in the inner regions of disks. In addition, this method may shed more light on the importance of radial drift in setting the inner disk H₂O reservoir by comparing small disks, where radial drift is expected to be very efficient, with large disks, in which drift is expected to be halted by substructures.

**Acknowledgements.** The authors would like to thank Andrea Banzatti for many very useful discussions. We also thank the referee for many thoughtful, constructive comments that helped improve the manuscript. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1282. The following National and International Funding Agencies funded and supported the MIRI development: NASA; ESA; Belgian Science Policy Office (BELSPO); Centre Nationale d'Etudes Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum für Luft- und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y Competitividad; Netherlands Research School for Astronomy (NOVA); Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Agency; and UK Space Agency. The authors acknowledge the GO #1640 team for developing their observing program with a zero-exclusive-access period, which includes CI Tau. M.T., E.v.D., and M.V. acknowledge support from the ERC grant 101019751 MOLDISK. E.v.D. acknowledges support the Danish National Research Foundation through the Center of Excellence “InterCat” (DNRF150). E.v.D., I.K., and A.M.A. acknowledge support from grant TOP-1 614.001.751 from the Dutch Research Council (NWO). D.G. thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA). B.T. is a Laureate of the Paris Region fellowship program, which is supported by the Ile-de-France Region and has received funding under the Horizon 2020 innovation framework program and Marie Sklodowska-Curie grant agreement No. 945298. T.H. and K.S. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. D.B. has been funded by Spanish MCIN/AEI/10.13039/501100011033 grants PID2019-107061GB-C61 and No. MDM-2017-0737. A.C.G. acknowledges from PRIN-MUR 2022 20228JPA3A “The path to star and planet formation in the JWST era (PATH)” funded by NextGeneration EU and by INAF-GoG 2022 “NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)” and Large Grant INAF 2022 “YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems (YODA)”. I.K. acknowledges funding from H2020-MSCA-ITN-2019, grant no. 860470 (CHAMELEON). G.P. gratefully acknowledges support from the Max Planck Society. This work also has made use of the following software packages that have not been mentioned in the main text: NumPy, SciPy, Astropy, Matplotlib, pandas, IPython, Jupyter (Harris et al. 2020; Virtanen et al. 2020; Astropy Collaboration 2013, 2018, 2022; Hunter 2007; Pandas Development Team 2020; Pérez & Granger 2007; Kluyver et al. 2016).

## References

Ardila, D. R., Basri, G., Walter, F. M., Valenti, J. A., & Johns-Krull, C. M. 2002, ApJ, 566, 1100

Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Avni, Y. 1976, ApJ, 210, 642

Banzatti, A., Pontoppidan, K. M., Salyk, C., et al. 2017, ApJ, 834, 152

Banzatti, A., Pascucci, I., Bosman, A. D., et al. 2020, ApJ, 903, 124

Banzatti, A., Abernathy, K. M., Brittain, S., et al. 2022, AJ, 163, 174

Banzatti, A., Pontoppidan, K. M., Carr, J. S., et al. 2023a, ApJ, 957, L22

Banzatti, A., Pontoppidan, K. M., Pére Chávez, J., et al. 2023b, AJ, 165, 72

Bast, J. E., Brown, J. M., Herczeg, G. J., van Dishoeck, E. F., & Pontoppidan, K. M. 2011, A&A, 527, A119

Bergin, E. A., & Williams, J. P. 2017, in Astrophysics and Space Science Library, 445, Formation, Evolution, and Dynamics of Young Solar Systems, eds. M. Pessah, & O. Gressel, 1

Blevins, S. M., Pontoppidan, K. M., Banzatti, A., et al. 2016, ApJ, 818, 22

Bosman, A. D., Tielens, A. G. G. M., & van Dishoeck, E. F. 2018, A&A, 611, A80

Bosman, A. D., Bergin, E. A., Calahan, J., & Duval, S. E. 2022, ApJ, 930, L26

Brown, J. M., Pontoppidan, K. M., van Dishoeck, E. F., et al. 2013, ApJ, 770, 94

Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2015, A&A, 575, A94

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, ASP Conf. Ser., 527, 583

Calahan, J. K., Bergin, E. A., & Bosman, A. D. 2022, ApJ, 934, L14

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Dawson, R. I., & Johnson, J. A. 2018, ARA&A, 56, 175

Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

Erb, D. 2022, https://doi.org/10.5281/zenodo.7255880

Facchini, S., van Dishoeck, E. F., Manara, C. F., et al. 2019, A&A, 626, L2

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gasman, D., van Dishoeck, E. F., Grant, S. L., et al. 2023, A&A, 679, A117

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, J. Quant. Spec. Radiat. Transf., 277, 107949

Grant, S. L., van Dishoeck, E. F., Tabone, B., et al. 2023, ApJ, 947, L6

GRAVITY Collaboration (Perraut, K., et al.) 2021, A&A, 655, A73

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Henning, T., Kamp, I., Samland, M., et al. 2024, PASP, 136, 054302

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Kalyaan, A., Pinilla, P., Krijt, S., Mulders, G. D., & Banzatti, A. 2021, ApJ, 921, 84

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt (IOS Press), 87

Lacy, J. H., Sneden, C., Kim, H., & Jaffe, D. T. 2017, ApJ, 838, 66

Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49

Meijerink, R., Pontoppidan, K. M., Blake, G. A., Poelman, D. R., & Dullemond, C. P. 2009, ApJ, 704, 1471

Morbidelli, A., Lunine, J. I., O'Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Annu. Rev. Earth Planet. Sci., 40, 251

Muñoz-Romero, C. E., Öberg, K. I., Banzatti, A., et al. 2024, ApJ, 964, 36

Mulders, G. D., Ciesla, F. J., Min, M., & Pascucci, I. 2015, ApJ, 807, 9

Najita, J. R., Carr, J. S., Strom, S. E., et al. 2010, ApJ, 712, 274

Najita, J. R., Carr, J. S., Salyk, C., et al. 2018, ApJ, 862, 122

Pandas Development Team 2020, https://doi.org/10.5281/zenodo.3509134

Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21

Perotti, G., Christiaens, V., Henning, T., et al. 2023, Nature, 620, 516

Pontoppidan, K. M., Salyk, C., Blake, G. A., et al. 2010, ApJ, 720, 887

Pontoppidan, K. M., Salyk, C., Banzatti, A., et al. 2024, ApJ, 963, 158

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in C. The art of scientific computing (Cambridge University Press)

Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

Salyk, C., Pontoppidan, K. M., Blake, G. A., et al. 2008, ApJ, 676, L49

Salyk, C., Pontoppidan, K. M., Blake, G. A., Najita, J. R., & Carr, J. S. 2011, ApJ, 731, 130

Salyk, C., Lacy, J., Richter, M., et al. 2019, ApJ, 874, 24

Schwarz, K. R., Henning, T., Christiaens, V., et al. 2024, ApJ, 962, 8

Tabone, B., van Hemert, M. C., van Dishoeck, E. F., & Black, J. H. 2021, A&A, 650, A192

Tabone, B., Bettoni, G., van Dishoeck, E. F., et al. 2023, Nat. Astron., 7, 805

Tabone, B., van Dishoeck, E. F., & Black, J. H. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202348487

Temmink, M., van Dishoeck, E. F., Grant, S. L., et al. 2024, A&A, 686, A117

Thi, W. F., van Dishoeck, E. F., Blake, G. A., et al. 2015, PASP, 127, 646

Wilson, T. L. 1999, Rep. Progr. Phys., 62, 143

Woitke, P., Min, M., Thi, W. F., et al. 2018, A&A, 618, A57

Wright, G. S., Wright, D., Goodson, G. B., et al. 2015, PASP, 127, 595

Zannese, M., Tabone, B., Habart, E., et al. 2024, Nat. Astron., 8, 577

Zhang, K., Pontoppidan, K. M., Salyk, C., & Blake, G. A. 2013, ApJ, 766, 82

Wells, M., Pel, J. W., Glasse, A., et al. 2015, PASP, 127, 646

Wilson, T. L. 1999, Rep. Progr. Phys., 62, 143

Woitke, P., Min, M., Thi, W. F., et al. 2018, A&A, 618, A57

Wright, G. S., Wright, D., Goodson, G. B., et al. 2015, PASP, 127, 595

Zannese, M., Tabone, B., Habart, E., et al. 2024, Nat. Astron., 8, 577

Zhang, K., Pontoppidan, K. M., Salyk, C., & Blake, G. A. 2013, ApJ, 766, 82
¹ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
e-mail: temmink@strw.leidenuniv.nl

² Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstraße 1, 85748 Garching, Germany

³ Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium

⁴ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France

⁵ Department of Astrophysics, University of Vienna, Türken-schanzstr. 17, 1180 Vienna, Austria

⁶ ETH Zürich, Institute for Particle Physics and Astrophysics, Wolfgang-Pauli-Str. 27, 8093 Zürich, Switzerland

⁷ Max-Planck-Institut für Astronomie (MPIA), Königstuhl 17, 69117 Heidelberg, Germany

⁸ Centro de Astrobiología (CAB), CSIC-INTA, ESAC Campus, Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

⁹ INAF – Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, 80131 Napoli, Italy

¹⁰ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, D02 XF86 Dublin, Ireland

¹¹ Kapteyn Astronomical Institute, Rijksuniversiteit Groningen, Postbus 800, 9700AV Groningen, The Netherlands

¹² Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
## Appendix A: Uncertainties for the JWST-MIRI subbands

<table><caption>Table A.1: Median continuum flux, estimated S/N from the ETC, and the obtained uncertainty σ for each subband.</caption><thead><tr><th>Subband ([µm])</th><th>Median flux [Jy]</th><th>ETC S/N</th><th>σ [mJy]</th></tr></thead><tbody><tr><td>1A (4.90-5.74)</td><td>1.31</td><td>533.5</td><td>2.5</td></tr><tr><td>1B (5.66-6.63)</td><td>1.31</td><td>613.5</td><td>2.1</td></tr><tr><td>1C (7.51-8.77)</td><td>1.32</td><td>743.9</td><td>1.8</td></tr><tr><td>2A (7.51-8.77)</td><td>1.43</td><td>745.2</td><td>1.9</td></tr><tr><td>2B (8.67-10.13)</td><td>1.98</td><td>962.6</td><td>2.1</td></tr><tr><td>2C (10.02-11.70)</td><td>2.03</td><td>1014.2</td><td>2.0</td></tr><tr><td>3A (11.55-13.47)</td><td>1.79</td><td>956.5</td><td>2.0</td></tr><tr><td>3B (13.34-15.57)</td><td>1.89</td><td>991.57</td><td>1.9</td></tr><tr><td>3C (15.41-17.98)</td><td>2.35</td><td>1143.0</td><td>2.1</td></tr><tr><td>4A (17.70-20.95)</td><td>2.69</td><td>497.3</td><td>5.4</td></tr><tr><td>4B (20.69-24.48)</td><td>2.85</td><td>271.7</td><td>10.5</td></tr><tr><td>4C (24.19-27.90)</td><td>2.97</td><td>75.5</td><td>39.4</td></tr></tbody></table>