© The Authors 2023

1 Introduction

The inner 0.1 to 10 au regions of protoplanetary disks are likely to be the cradle for terrestrial planets around low-mass stars. The high temperatures (≥100 K) and densities (≥10⁸ cm⁻³) in these regions, along with the locations of the H₂O and CO₂ snow lines and the presence of substructures (see e.g. Grant et al. 2023; Tabone et al. 2023; van Dishoeck et al. 2023; Pontoppidan et al. 2014 for a review), dictate the composition of the gas and therefore the elemental abundances available to atmosphere formation of accreting planets.

The mid-infrared wavelength range observed by the Medium Resolution Spectroscopy (MRS; Wells et al. 2015; Argyriou et al. 2023) mode of the Mid-InfraRed Instrument (MIRI; Wright et al. 2015, 2023; Rieke et al. 2015) on board the James Webb Space Telescope (JWST; Rigby et al. 2023) allows us to examine these inner regions of protoplanetary disks. Its wide wavelength range (4.9 to 28.1 µm) covers a variety of molecular features, including the large forest of H₂O lines: from the ro-vibrational bending mode from 5 to 8 µm, to the pure rotational lines around 10 µm and onwards (Meijerink et al. 2009). These lines are thought to probe radially different regions of the disk. Generally, the temperature of the gas probed decreases as the wavelength increases, likely corresponding to moving from the inner disk outwards (e.g. Banzatti et al. 2017, 2023). The inner disk chemistry is now starting to be seen with MIRI/MRS (e.g. Kóspál et al. 2023; Grant et al. 2023; Tabone et al. 2023; Kamp et al. 2023; van Dishoeck et al. 2023; Perotti et al. 2023).

Excitation of H₂O can occur due to collisions with H, H₂, He, and electrons; radiation from hot dust; photodesorption from dust grains; and chemical formation (Meijerink et al. 2009; van Dishoeck et al. 2013, 2021). Already in the era of Spitzer, its InfraRed Spectrograph (IRS) unveiled a large diversity between the compositions of disks around T Tauri stars (e.g. Pontoppidan et al. 2010, 2014; Carr & Najita 2011). The samples examined by Carr & Najita (2011) and Najita et al. (2013) showed that the ratio of HCN to H₂O increases with increasing disk mass. Recently, Banzatti et al. (2020) observed a similar correlation: the ratio of H₂O versus carbon-bearing molecules seems to correlate with disk size. This is interpreted as a smaller disk size indicating more efficient drift of icy pebbles, allowing the inner disk to be replenished with an ice reservoir that may sublimate. On the other hand, substructures in larger disks can prevent this transport. The disk of Sz 98, which is studied here, contains both of these features, being a large dust disk with rings, but also showing bright millimetre wavelength emission within several tens of au from the star, which we refer to as the central core (e.g. van Terwisga et al. 2018).

A pathway for formation of H₂O in the gas phase is from its precursor OH: OH + H₂ → H₂O + H, which is most efficient at higher temperatures, typically >250 K (e.g. Woitke et al. 2009; Glassgold et al. 2009; van Dishoeck et al. 2013). In reverse, the OH reservoir can be replenished again by photodissociation of H₂O, which puts OH back into the cycle (e.g. Harich et al. 2000; van Harrevelt & van Hemert 2000; Tabone et al. 2021). The balance of the pebble drift and chemical formation from OH producing H₂O, versus destruction by irradiation of the disk ultimately dictate the abundance of gaseous H₂O across the disk probed here in the infrared.

Sz 98 is a relatively cool, actively accreting classical T Tauri star (Merín et al. 2008; Mortier et al. 2011) of spectral type K7 (Alcalá et al. 2017) with a disk amongst the 2% largest and brightest dust disks in Lupus (van Terwisga et al. 2018). It has a luminosity L_* = 1.5 L_⊙ and a mass M_* = 0.74 M_⊙ (Alcalá et al. 2017) at a Gaia distance of approximately 156 pc (Gaia Collaboration 2016, 2023), and disk radii R_gas and R_dust of 360 and 180 au, respectively (Ansdell et al. 2018). The mass of the disk is estimated to be 0.07 M_⊙ (van Terwisga et al. 2019). The star has an accretion rate and luminosity of Ṁ_acc ≈ 10^−7.6 M_⊙yr⁻¹ (Alcalá et al. 2017) and L_acc ≈ 10^−0.72 L_⊙ (Nisini et al. 2018). In general the disk has not been found to have large substructures, aside from a small continuum break around 80 au and a ring around 90 au (Tazzari et al. 2017; van Terwisga et al. 2018; van der Marel et al. 2019; Miotello et al. 2019). Some evidence for additional ring-like structure has been found around ~120 au (van Terwisga et al. 2019). Furthermore, based on scattered light images from the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Garufi et al. 2022) instrument, the disk’s inner rim might be casting a uniform shadow on the outer disk. Miotello et al. (2019) suggest that the outer disk is depleted in volatile gaseous carbon and oxygen, and find signs for a gaseous C/O ratio larger than unity based on bright C₂H and faint ¹³CO emission seen in Atacama Large Millimeter/submillimeter Array (ALMA) observations. The dust grains are expected to have grown at least up to millimetre sizes (Lommen et al. 2007; Ubach et al. 2012). The dust within 0.5 au from the star is likely warm (dust excess of T ~ 600 K) based on photometry with Spitzer (Wahhaj et al. 2010). Sz 98 was previously observed with Spitzer/IRS in 2008, but in its low resolution mode only. The higher sensitivity and spectral resolving resolution of the MRS ranging from a resolving power λ/Δλ ≈ 3500 at shorter wavelengths to λ/Δλ ≈ 1500 at the longer wavelengths (Jones et al. 2023) now allows us to detect H₂O in all of its bands.

Since we detect water across the full MIRI/MRS wavelength range, the focus of this paper is on the H₂O. The paper is structured as follows. The data reduction is described in Sect. 2, where we also discuss our modelling methods. The resulting spectrum and best model fits are shown in Sect. 3. A discussion regarding the implications on the chemistry in the inner disk of Sz 98 and a comparison to other disks can be found in Sect. 4. Finally, we summarise our conclusions in Sect. 5.

2 Methods

2.1 Data acquisition and reduction

Sz 98 was observed with MIRI/MRS as part of the MIRI mid-INfrared Disk Survey (MINDS) JWST GTO Programme (PID: 1282, PI: T. Henning). Taken on August 8 2022, the exposure of all twelve bands in a four-point dither pattern resulted in an exposure time of approximately 14 min per grating setting.

The data were processed using version 1.9.4¹ of the JWST pipeline (Bushouse et al. 2022). We followed the methods described in Gasman et al. (2023) that are specific to point sources to defringe and flux calibrate the data. The reference files exploit the repeatability of the fringes when the pointing is consistent between the science observation and the reference target, and are extracted from the A star HD 163466 (PID: 1050). In the Sz 98 case a more significant pointing offset is seen in channel 2 (~7.5–11.7 µm), which is dominated by the silicate feature, but we do not analyse this region in detail. Similarly, a clean spectrophotometric calibration is derived on-sky from HD 163466. This method is adopted mainly for the advantages related to defringing, since the standard pipeline fringe flats require the additional residual_fringe step to defringe the spectrum. As shown by Gasman et al. (2023), the latter step may change the shape of molecular features, which would affect the derived excitation properties.

In order to extract the spectrum, the signal in an aperture of 2.5 × FWHM centred on the source was summed. The background was estimated from an annulus that linearly grows between 5 × FWHM and 7.5 × FWHM in the shortest wavelengths, and 3 × FWHM to 3.75 × FWHM in the longest wavelengths. The aperture correction factors applied are the same as those in Argyriou et al. (2023), and account for the fraction of the signal of the Point Spread Function (PSF) outside of the aperture, and inside the annulus. Additionally, the outlier detection step of spec3 was skipped, since this has spurious results for the data taken with the MRS due to undersampling of the PSF. No rescaling of sub-bands was required to stitch the spectrum.

2.2 Slab models and fitting procedure

We subtracted the continuum from the spectrum by fitting a cubic spline through the line-free sections, in order to fit local thermodynamic equilibrium (LTE) slab models. Line-free sections were selected iteratively from visual inspection of the slab fits. Especially in the region of 24 µm and beyond, the data become noisy and filled with artefacts caused by the low signal in the reference A star HD 163466. Using this knowledge, the continuum in the >24 µm region was selected in such a way that artefacts are avoided. To do so, the continuum points were placed to make the continuum follow along with the larger artefacts where the signal of the A star drops close to 0 (see Gasman et al. 2023 for more details on how the reference files were derived). Due to the presence of these artefacts, the line fluxes in this region may have a larger uncertainty. However, since we do not analyse this region in great detail due to fewer features of interest, this does not affect our results.

For the H₂O lines longwards of ~10 µm, LTE is a good first estimate, though this assumption is less accurate for the lines around 6.5 µm (e.g. Meijerink et al. 2009; Bosman et al. 2022a). The strongest deviation from LTE will be seen in the high-energy lines (E_up > 3000 K, e.g. Meijerink et al. 2009). The line profiles were assumed to be Gaussian, with a broadening of ΔV = 4.7 km s⁻¹ (σ = 2 km s⁻¹), similarly to Salyk et al. (2011). For molecules with densely packed lines, in other words the CO₂ and HCN Q-branches, we included mutual shielding from adjacent lines as described in Tabone et al. (2023). Subsequently, three parameters were varied in order to fit the features in the observation: the column density N, the excitation temperature T, and the emitting area πR². The latter scales the strength of the features to match the strength in the observed spectrum. We note that the excitation temperature does not need to be the same as the kinetic temperature of the gas.

The slab models were convolved to a constant resolution per region in the same order as that of the MRS in the relevant bands (ranging from a resolving power λ/Δλ ≈ 3500 at shorter wavelengths to λ/Δλ ≈ 1500 at the longer wavelengths, specific values given in Table 1, Jones et al. 2023), and resampled with spectres (Carnall 2017) to the same wavelength grid. Since we examined the H₂O features over a large range of wavelengths, the resolution greatly varied between different fits. Additionally, the regions of the disk probed by the different wavelengths and their excitation conditions are not expected to be constant throughout, hence we treated the shorter wavelengths separately from the longer wavelengths (e.g. Banzatti et al. 2023). By dividing the spectrum into the 5–6.5 µm region, the ~13.6–16.3 µm region, the ~ 17–23 µm region, and 23 µm and onwards, regions of similar excitation conditions were addressed separately.

Using χ² fitting, the most likely values of the variables are found. Similarly to Grant et al. (2023), we first identified relatively bright and isolated lines, used these to find the best fitting slab model, subtracted the model, and moved on to other molecular features. Due to the spectra being very dense in lines, this became an iterative procedure where the noise estimates were taken from the line-subtracted spectra. The order in which we fit the molecules per region is as follows: region 1 – CO, H₂O; region 2 – H₂O, OH, HCN, CO₂; region 3 – H₂O, OH; region 4– H₂O, OH. The emitting area is parametrised in terms of an arbitrary disk radius R. We note that this does not need to be the radius at which the emission is located, but rather the equivalent emitting area. For example, the emission could be confined within a ring of a total emitting area of πR². The number of molecules 𝒩 is also included, resulting from the column density and emitting area. It is a more robust metric for optically thin species.

Similarly to Grant et al. (2023), the reduced χ² is then defined using: $${\chi }^2=\frac{1}{M}{\displaystyle \sum _{i=1}^M\frac{F_{\text{obs,}i}-F_{\text{mod.}i}}{{\sigma }^2}},$$(1)

where M is the number of data points in the selected wavelength window, σ the noise estimated from the selected region with the lines removed, and F_obs and F_mod are the observed and modelled continuum-subtracted flux, respectively. The confidence intervals are defined as ${\chi }_{\min }^2+2.3$, ${\chi }_{\min }^2+6.2$, and ${\chi }_{\min }^2+11.8$; for 1σ, 2σ, and 3σ, respectively (Avni 1976; Press et al. 1992). The σ per region was estimated from a standard deviation on the spectrum itself, after subtracting the best-fit slab models. This was done since, although the visually most ‘line-free’ regions were selected, faint lines were often still present. The resulting spectra and the noise regions are given in Sect. 3, along with the regions used to fit the lines.

The spacing of the grid is consistent between molecules (ΔT = 25 K, Δ log₁₀(N) = 0.16 cm⁻², and Δ log₁₀ (R) = 0.02 au), but the range was varied depending on how hot or cold the excitation is expected to be. For OH much higher temperatures (500–4000 K) were used compared to CO₂ (100–1400 K) and H₂O (100–1500 K).

Fig. 1 Full MIRI/MRS spectrum of the Sz 98 disk. The extent per region is indicated, along with the detected species. The Spitzer LRS spectrum from the CASSIS database is shown in light grey (Lebouteiller et al. 2011).

3 Results

We present the full spectrum in Fig. 1, where we indicate the molecules detected per region. The overall shape of the continuum is typical of a T Tauri disk: a discernible silicate feature around 10 µm indicating the presence of small silicate grains, and excess emission at the longer wavelengths. The silicate feature around 10 µm is sensitive to changes in grain size, where larger grain sizes indicate more evolved dust with a lower opacity (e.g. Bouwman et al. 2001; Przygodda et al. 2003; Kessler-Silacci et al. 2006; Juhász et al. 2010). The peak value normalised to the continuum is ~1.9, which indicates a relatively small grain size of a few µm based on the models in Kessler-Silacci et al. (2006). The peak value seems to be slightly below the peak value of EX Lup, which has since been observed with MIRI/MRS as well (Kóspál et al. 2023).

On top of the continuum, we detect CO around 5 µm, CO₂, HCN, OH from ~13 µm and onwards; and, most strikingly, we see H₂O features from 5 µm up to 27 µm. The relative strength of emission for CO₂ and H₂O is opposite to the GW Lup case (Grant et al. 2023), where CO₂ was found to be much stronger than H₂O. Additionally, there is a general lack of detectable carbon-bearing species in the inner disk: aside from CO₂ and HCN, we detect no organic molecules. Most notably, C₂H₂ is not detected.

The best-fit parameters are presented in Table 1, and the corresponding fits per region overlaid on the data can be found in Fig. 2. The σ per region are documented in both Table 1 and Fig. 2. The χ² maps representing the confidence per fit are included in Appendix A. Some sections of the continuum subtracted spectrum are negative, particularly in spectral region 4. This is due to the selection of the continuum, which is greatly influenced by the noise and artefacts in this part of the spectrum. We note an up to ~30% flux discrepancy between the Spitzer/IRS and MRS continuum shortwards of ~16 µm. Alcalá et al. (2017) found a similar discrepancy between photometry data and X-shooter spectroscopy (taken in 2015) of the object, and noted that this was within the expected variability range for Class II young stellar objects found by Venuti et al. (2014), Fischer et al. (2023). Additionally, Bredall et al. (2020) classified Sz 98 as a ‘dipper’ star, which typically thought to be caused by the disk being close to edge-on (e.g. Stauffer et al. 2015; Bodman et al. 2017). However, based on ALMA millimetre emission, its inclination is 47.1° (see Tazzari et al. 2017 and Appendix D). It is therefore possible that the inner and outer disk are misaligned, causing the discrepancy to be different in the shorter and longer wavelengths.

Despite being located in similar spectral regions, a wide variety of best fit parameters is found for different species. Due to this, we note that it is unlikely they all emit from the same region, as some species may be located deeper or farther out in the disk, and different energy levels are probed per species. Furthermore, the actual spectral overlap is minimal, despite what the MRS resolution might imply. A simple addition of the different models per region is therefore a good approximation. Additionally, simple LTE excitation may not be the most fitting assumption for some species. In the following sections we discuss the bestfit results for the detected molecules in more detail, and some non-detections.

3.1 H₂O

In Fig. 2, we present sections of the continuum-subtracted spectrum where pronounced H₂O features are present, with the best-fit results of different species. As noted in several previous works (e.g. Blevins et al. 2016; Banzatti et al. 2023), it is expected that the inner (warmer) to outer (colder) regions of the inner disk are probed from shorter to longer wavelengths, respectively. Indeed, we can conclude that this is the case for the spectrum of Sz 98, based on the best-fit parameters presented in Table 1, and the confidence intervals shown in Fig. A.1. The best-fit temperature and emitting radius change as we move to longer wavelengths. The temperatures slowly decrease from 950 to 250 K, indicating that we are gradually probing colder and/or less excited gas. We present this tentative trend in Fig. 3. The error bars are based on the 1σ-contours of the χ² plots in Fig. A.1. The best-fit slab model of region 1 will underestimate the lines in region 4, and vice versa. Furthermore, adding the H₂O spectra from all regions together significantly overestimates the flux. In reality a specific disk region of a certain temperature and radial extent will not be contained to a specific spectral region, but influences lines in other regions as well. Alternatively, adding the spectra together by assuming the radii found are instead the inner and outer radii of a series of annuli, the spectrum is better reproduced while keeping the radii within the 1σ-confidence intervals. This indicates that one slab cannot be used to fit the entire MIRI/MRS range, and future work must accommodate for a temperature and emitting area gradient in the fitting procedure.

Additionally, we systematically find the H₂O lines to be optically thick, and our H₂O column densities are in a similar range as those found previously from Spitzer spectra of other disks (Carr & Najita 2011; Salyk et al. 2011), although in regions 3 and 4 they are higher than usually inferred. As demonstrated by the models of Meijerink et al. (2009); Walsh et al. (2015), we are likely not probing the full column density of H₂O down to the mid-plane, especially in the short wavelength region. Part of the H₂O gas is hidden below the dust continuum where τ_IR ≈ 1, where dust might be blocking emission from the deeper layers in the disk. The small emitting radius of H₂O in spectral region 1 indicates that we would be probing the inner gas disk (e.g. Dullemond & Monnier 2010).

Not all H₂O lines are fit equally well. Meijerink et al. (2009) show that high energy lines are less likely to be thermally excited, and might be better fit with lower temperature slab models. This would be a sign that some lines are sub-thermally excited, and the LTE assumption is not applicable. However, the LTE slab models fit the spectrum very well, therefore no significant evidence for non-LTE excitation of H₂O is found here. While more detailed thermo-chemical models could result in a better representation of the spectrum, the increased complexity introduces more uncertainties in the fits, and this is left for future work, along with fits of temperature gradients. However, the temperature trend seen in the simple slab models of the different regions is robust.

In order to assess the correctness of the χ² fit results, we compare the fluxes of pairs of lines with the same upper energy levels, but different Einstein A_ul coefficients. The flux ratio of these lines will depend primarily on the opacity of the lines rather than the temperature, providing a robust estimate for the column density. To reduce the effects of artefacts and other molecular features, we limit the examined range to isolated lines. In this manner, two pairs of lines were identified that are primarily sensitive to the changes in column density. The properties of these lines can be found in Table 2. Evaluating the ratios of these lines in slab models of changing column density and temperature results in the coloured trends in panels a and b of Fig. 4, while the black horizontal lines result from the flux ratios in the data. For low column densities below ~10¹⁶ cm⁻²−10¹⁸ cm⁻² (depending on the temperature assumed), the trends are largely flat, since the flux ratio depends on the A_ul ratio. Once one of the lines becomes opaque, this ratio will change, resulting in the upward trend for higher column densities. When assessing the observed flux ratios in the data the column density should indeed be high as suggested by the best-fit slab models. For optically thin lines, the flux ratios converge to a single value that corresponds to the A_ul ratio. This is not the case, as shown in Table 2, indicating that the brighter line could be opaque. However, we note that these lines are part of a cluster of lines, and the ratios are likely affected by line blending. In that case, we cannot claim that the lines are optically thick based on the fact that the ratio of A_ul does not equal the flux ratio. However, the flux ratios of the slab models presented in Fig. 4 are similarly affected, therefore these trends are still representative.

Note that it is assumed that the emitting area of the line pairs is equal. Assuming the longer wavelength line traces a larger area, the larger the discrepancy between the emitting areas, the higher the flux of the reference line at longer wavelengths, and the more the trends are shifted down. Therefore, even for differing emitting areas, this conclusion remains valid.

Since H₂¹⁸O-rich regions are expected to be located deeper in the disk atmosphere emitting optically thin lines, detection and characterisation of the isotopologue would put constraints on the H₂O column densities, and allow for the abundance to be measured (Calahan et al. 2022). Several isolated lines are potentially detectable by MIRI/MRS, which we plot in Fig. 5. Given the large column density of H₂O, the isolated H₂¹⁸O could potentially be visible. We compare the data to slabs of the same excitation temperature and emitting radius as the bestfit H₂O values in the wavelength region assuming a ¹⁶O/¹⁸O of 550. In Fig. 5, we demonstrate that H₂¹⁸O is not detected in our data. One of the brightest, isolated lines in the spectrum is the line around 27 µm, the details of which can be found in Table 2. Since it is not detected, the lower limit on the line flux of the H₂¹⁸O line (taken from 26.9626 to 27.0106 µm) is ~0.4×10⁻¹⁴ erg s⁻¹ cm⁻² normalised to a distance of 140 pc. Assuming the integrated flux of the H₂¹⁸O feature is in this upper limit on the line flux, the maximum column density of H₂O is 3.5 × 10¹⁹ cm⁻² when assuming the ISM H₂¹⁶O/H₂¹⁸O ratio of 550. This value has been indicated as an upper limit in Fig. 4, and in the χ² maps in Appendix A. Based on the non-detection, we can plot a comparison to the slab model line ratios in a similar fashion as for H₂O. However, now the comparison is between the H₂¹⁸O line and a H₂¹⁶O line of similar upper energy (see Table 2 for the details), shown in panel c of Fig. 4. Note that the column density on the x-axis and in the slab model of H₂¹⁶O is multiplied by 550 compared to the column density of H₂¹⁸O. Figure 4 indicates that the column densities are indeed high, but likely between 5 × 10¹⁸ to 3.5 × 10¹⁹ cm⁻² rather than 7.5 × 10¹⁹ cm⁻² as found for region 3 in Table 1. The larger column density might be possible, provided the H₂¹⁶O/H₂¹⁸O ratio is larger than 550. Some variation in this ratio per transition might be possible (see e.g. Calahan et al. 2022).

Fig. 2 Slab fits for the four different wavelength ranges. The top and bottom panels per region are the data and total model, and the individual slab models; respectively. Spurious spikes and features from the data reduction have been blanked out (see the list in Appendix B). The horizontal lines indicate the regions used to fit the molecules, or the region where σ is estimated after subtracting the slab models.

Fig. 3 Best-fit temperature and emitting radius for H2O per region of the inner disk. The error bars indicate the 1σ contours of the χ2 maps.

Fig. 4 Flux ratios of H216O lines of the same upper level in the data compared to slab models (top); and H216O/H218O ratio in the data and slab models (bottom). The properties of the transitions can be found in Table 2. The dotted black vertical line indicates the upper limit of the column density based on the 1σ level of the spectrum and the nondetection of the H218O line. The x-axis shows the column density of H216O. In the bottom panel, H216O/H218O ratio of 550 is assumed.

3.2 CO

The best-fit parameters for the CO emission in the disk are 1675 K, with a radius of 1.00 au, and a column density of 1.4× 10¹⁵ cm⁻². However, we note a degeneracy between the three parameters, causing the confidence intervals to be relatively wide. The best fit excitation temperature found here is not likely equal to the true gas temperature, which is likely colder than the best fit of 1675 K and optically thick, since the gas temperatures probed by CO are typically not more than 1000 K (e.g. Brown et al. 2013; Banzatti et al. 2022; Anderson et al. 2021). As shown in Table 1, the upper energy levels in this wavelength region are typically much higher than that of H₂O; therefore, we are likely not tracing the same region for both species. Only the rotationally excited CO lines of the υ = 1–0 band can be detected by MIRI/MRS. The CO lines in the Near Infrared Spectrograph (NIRSpec) region and observed from the ground for bright sources can be expected to be probing more similar conditions.

3.3 CO₂

We detect the CO₂ Q-branch, although the hot bands are not confidently detected, similarly to the isotopologue ¹³CO₂. The CO₂ emission is much weaker compared to the H₂O emission, which is opposite to the case of GW Lup (Grant et al. 2023). Due to the blending with H₂O lines and faintness of the features, it is difficult to constrain the excitation properties. A demonstration of the differences between optically thick and optically thin emission with ¹³CO₂ can be found in Appendix C. The best-fit temperature when including the hot band around 13.8 µm is relatively cold, with 125 K from an emitting radius of 1.63 au, and a column density of ~2.4× 10¹⁹ cm⁻². Due to the lower excitation temperature, it likely originates from a deeper layer in the disk, or farther away from the star. However, when fitting only the Q-branch, the best fit parameters are an excitation temperature of 250 K, and emitting radius of 0.28 au (see Fig. A.2). While the emitting area in this case is similar to that of H₂O, the excitation temperature is different.

3.4 OH

The maximum E_u of the OH transitions we tentatively detect is ~15 000 K, corresponding to an upper rotational quantum number of ~25 (Tabone et al. 2021). The OH emission, though poorly constrained, seems to have a high excitation temperature (likely nearer 2000 K or higher) with a lower column density than H₂O (<2 × 10¹⁷ cm⁻²) in all regions where it is detectable. Under the LTE assumption, the emission could originate from high in the disk atmosphere, based on these elevated temperatures. However, other explanations are possible. We further discuss the excitation of OH in this context in Sect. 4.5.

Fig. 5 Different regions where H218O is expected to show lines sufficiently isolated from H2O. These lines are not seen in the spectrum, compared a H218O slab models assuming the same best-fit parameters as H2O in these regions, assuming a 16O/18O of 550. The models have been multiplied by 5 for visibility. The H2O slab models are based on the best-fit parameters presented in Table 1 per wavelength region.

3.5 HCN

The broad HCN Q-branch is detected in spectral region 2, around 14 µm; along with a hot band around 14.3 µm. The constraints on the fit are poor due to the degeneracy between the column density and the emitting radius for low column density, and the degeneracy between the column density and temperature for high column density. While HCN can be fit with both an optically thin and optically thick solution, we find a best fit for a temperature of 1075 K, with a radius of 0.76 au, and a column density of 4.3 × 10¹⁶ cm⁻². However, based on thermo-chemical models, it is more likely in the range of 330 K at a much higher column density (Woitke et al. 2018), which is still contained in the 1σ-confidence level (see Fig. A.2). This is far higher than previously found from Spitzer data (e.g. Salyk et al. 2011; Carr & Najita 2011), which could be an incorrect fit resulting from the lower resolution. For example, Grant et al. (2023) showed that the column density of CO₂ inferred from MIRI/MRS data must be much higher than previously assumed for GW Lup.

3.6 Hydrogen

We report the detection of a hydrogen recombination line. The strongest hydrogen recombination line, HI (6−5) at ~7.5 µm, is notably weaker than in GW Lup, where [NeII] was also detected (Grant et al. 2023). The location of this line is indicated in Fig. 1. Surprisingly, the HI (7−6) line, which is thought to trace accretion, is not detected or hidden due to a H₂O line (Rigliaco et al. 2015), despite Sz 98 being an active accretor (Merín et al. 2008; Mortier et al. 2011).

3.7 Non-detections

Notably, some species are not detected currently. Among these are the hydrocarbons, C₂H₂ and CH₄, of which the former has been detected in GW Lup (Grant et al. 2023). NH₃ is also not present. Despite being a strong accretor, no [NeII] is found. Furthermore, no molecular hydrogen is found to be strong enough to be detected in the forest of H₂O lines. As mentioned above, while H₂O is easily visible in the spectrum, its isotopologue H₂¹⁸O is not.

4 Discussion

As presented in Sect. 3, one of the defining features of the spectrum of Sz 98 is that H₂O emission is dominant, compared to relatively weak CO₂ and a lack of C₂H₂, together indicating a C/O ratio below unity. This is surprising, considering the large disk size and presence of dust traps, which should indicate limited drift of icy pebbles and thus correlate with a relative increase in carbon-bearing molecules (Banzatti et al. 2020). However, this correlation depends on assumptions related to the ice composition, formation timescale of substructures, among others; and may therefore change with time (e.g. Piso et al. 2015). In contrast, Miotello et al. (2019) inferred that the volatile C/O ratio is >1 in the outer disk of Sz 98 based on bright C₂H emission, inconsistent with our findings for the gas composition of the inner disk. We discuss a number of processes that may explain the observation of the inner disk of Sz 98.

Fig. 6 ALMA millimetre dust continuum emission from Sz 98 centred at 232.984 GHz (1286.75 µm). The dotted lines indicate gaps (G), the solid lines indicate rings (R), and their radii are included in au. Appendix D describes how the image is created.

4.1 Radial drift increasing the H₂O-gas reservoir

Despite being a very large and extended dust disk, the millimetre continuum shows a considerable central concentration resulting in the characteristic ‘fried egg’ shape found in van Terwisga et al. (2018) and presented in Fig. 6. While various explanations are possible for this characteristic shape, it could be indicative of radial transport inwards. If the bright central continuum is indeed due to pebble drift, this radial transport means pebbles containing H₂O-ice fall inwards within the H₂O snow line, estimated to lie at ~1.1 au (see Appendix D). This results in an increase in the volatile oxygen reservoir due to sublimation of H₂O ice (e.g. Ciesla & Cuzzi 2006; Banzatti et al. 2020). The inner disk is expected to be more replenished in oxygen-bearing species compared to carbon due to the snow lines of CO and CO₂ being farther out than that of H₂O (e.g. Öberg & Bergin 2021). The CO₂ and CO snow lines are estimated to lie at 2.2 au and 20 au, respectively, in Appendix D. Therefore, despite not conforming with the disk mass or size versus carbon bearing species correlation (Carr & Najita 2011; Najita et al. 2013; Banzatti et al. 2020), radial transport could still be the cause of the large H₂O column in the inner disk. This could indicate that the substructures of the disk are either not 100% efficient at trapping pebbles (see e.g. Sturm et al. 2022), or were formed late in the disk evolution when the pebbles containing H₂O-ice had already drifted inwards. The age of Sz 98 has been estimated to be between 1.7 and 5.6 Myr (van der Marel et al. 2019), and typical drifting timescales can be anywhere between ~ 10 kyr and ~ 1 Myr starting shortly after disk formation (Birnstiel et al. 2012, 2015). On the other hand, the timescale on which substructures form is uncertain, although relatively young objects exist that already show gaps, such as HL Tau (ALMA Partnership et al. 2015) (~1 Myr old, van der Marel et al. 2019). Therefore, pebbles could have drifted inwards prior to the formation of the gaps, but more information about gap formation is required to confirm this.

4.2 H₂O across the inner disk

The H₂O lines in the mid-infrared are thought to probe the ~0.1 au region out to the H₂O snow line (e.g. Banzatti et al. 2017, Banzatti et al. 2023). The critical densities to excite H₂O emission are higher (~10¹³ cm⁻³) at shorter wavelengths for the vibrational bending modes; and lower (~10⁸ cm⁻³) for the rotational lines at the longer wavelengths (Meijerink et al. 2009). This indicates that excitation due to collisions is less efficient for the vibrational bending modes, requiring a larger density for thermal excitation. Similarly, the rotational lines are more easily excited, requiring lower density. Indeed, in the shorter wavelengths around 3–4 µm, previous studies using a variety of ground-based instruments have found high temperature (1500 K), optically thick (10²⁰ cm⁻²) H₂O emission in other disks (Carr et al. 2004; Salyk et al. 2009, 2022; Doppmann et al. 2011; Banzatti et al. 2023). Moving outwards to the longer wavelengths probed by the MRS, the emission is well fit with 300–600 K with a column density of 10¹⁸ cm⁻² based on Spitzer IRS (e.g. Carr & Najita 2011; Salyk et al. 2011).

These ranges are not unlike what we find here (see Table 1): the ~5 µm range shows hotter emission around 950 K, while the best-fit temperatures slowly decrease down to 250–650 K at the longer wavelengths; although we find higher column densities in regions 3 and 4. Additionally, the column density also shows the opposite trend, being higher for regions 2 and 3 compared to region 1. However, we note that the maps in Fig. A.1 show that the fits are less well-constrained in terms of column density, and could all very well be in the 10¹⁸ cm⁻² range still. Alternatively, we could be probing a change in dust opacity (Antonellini et al. 2015), where τ_dust is located higher up in the disk closer to the star. Additionally, Banzatti et al. (2023) find the slab model fits to result in similar column densities for different sources throughout different wavelengths, of the order of ~10¹⁸ cm⁻². They therefore suggest that the different wavelengths probe the conditions where excitation is met per disk radius, resulting in decreasing temperatures and increasing emitting radii, at similar column densities. Since the column density is not well-constrained, this may still be true here.

Finally, some disks in the sample of Pontoppidan et al. (2010) exhibit relatively lower line fluxes at longer wavelengths in Spitzer, pointing to some depletion past the mid-plane snow line, due to the ‘cold finger effect’, which results in vertical transport across the snow line. In Appendix D, the snow line is found to be located at 1.1 au. Due to the brightness of the lines in region 4, it is unlikely that a significant amount of H₂O is transported down to the mid-plane across the H₂O snow line in Sz 98, since this would result in relatively weaker lines in region 4 due to formation of H₂O-ice (see also Pontoppidan et al. 2010). This indicates that a high abundance of H₂O is present across the disk surface. Furthermore, if the emitting radius of 1.4 au for region 4 found in Table 1 corresponds to the true emitting radius, it is indeed close to the 1.1 au snow line. However, far-IR H₂O past the MRS wavelength range would be more telling regarding the presence or absence of a ‘cold finger effect’, as done in Blevins et al. (2016).

4.3 H₂O self-shielding

H₂O is capable of self-shielding against UV radiation when it has sufficiently high column densities, starting at ~2 × 10¹⁷ cm⁻² (Bethell & Bergin 2009; Heays et al. 2017). The column densities inferred here (see Table 1), are well above this range. In this case, UV radiation cannot penetrate deep into the H₂O column, resulting in less photodissociation of H₂O, and potentially shielding other species (Bosman et al. 2022a).

As mentioned above, OH can be recycled back from H₂O through photodissociation of H₂O. Therefore, H₂O self-shielding results in a decreased OH abundance. The upper disk will be exposed to irradiation, resulting in a relatively thin OH layer above the self-shielding H₂O column (Walsh et al. 2015). Therefore, the gas containing OH is expected to be hot, with a lower column density. Assuming the excitation is in LTE, the fact that high energy lines are clearly observed at least in the ~ 13–16 µm range indicates that the OH gas must be hot. Indeed, from our fit (see Table 1) this is what we find. However, both prompt emission and chemical pumping from O + H₂ can also excite these higher energy lines. We discuss this further in Sect. 4.5

4.4 H₂O and OH line fluxes

The OH emission is relatively weak in Fig. 2 compared to H₂O. Banzatti et al. (2017) examine the line fluxes of the 12.52 µm and the tentatively detected 12.6 µm H₂O and OH lines, respectively, of a sample of disks observed with Spitzer. We plot the same line fluxes from these Spitzer/IRS spectra compared to the accretion luminosity in Fig. 7, where in the top panel only the H₂O line flux is shown, and the bottom panel the OH/H₂O-ratio. For samples where only an upper limit can be given for both H₂O and OH, the ratio is not included. The line fluxes for Sz 98 are calculated from the MIRI/MRS spectrum over the same interval, with the best-fit slab models of other species subtracted. Generally, the line flux of H₂O increases with stellar luminosity (Salyk et al. 2011), and stellar mass. From the top panel of Fig. 7, it becomes clear that the flux of the 12.52 µm line is on the lower end for Sz 98, but not out of the ordinary. A lower line flux may indicate that the inner disk has a comparatively higher amount of small dust blocking part of the H₂O column (Antonellini et al. 2017). In this case part of the radiation can be extincted by the dust, rather than the H₂O itself, making H₂O self-shielding less important for the chemistry in the inner disk. Additionally, the line flux ratio with OH shows a similar trend. The line flux of OH compared to H₂O is relatively low, but not unique in the sample.

Fig. 7 Line fluxes of H2O (top) and line flux ratio OH/H2O (bottom) and accretion luminosities from the disk sample in Banzatti et al. (2017) compared to Sz 98 (red dot). The size of the circle is proportional to the stellar mass, and the arrows indicate upper and lower limits. All fluxes have been normalised to 140 pc.

4.5 OH prompt emission

H₂O photodissociation by UV photons with a wavelength of <144 nm is known to produce OH in high rotationally excited states (Harich et al. 2000; van Harrevelt & van Hemert 2000). Tabone et al. (2021) show that newborn OH formed by H₂O photodissociation produces a series of rotationally excited lines detectable longwards of 9 µm, a process called prompt emission. In the MIRI/MRS, this process can be traced by highly excited lines shortwards of 10 µm that can uniquely be excited by H₂O photodissociation. Tabone et al. (2021) show that the series of rotational lines of OH exhibit relatively constant photon flux. Based on the strength of our detected OH lines longwards of 14 µm, the lines in the 10 µm would be detected, but already around 12.65 µm the lines are much weaker. Figure 8, shows that only very weak, if any, OH emission is found in this wavelength region. In case of prompt emission, it is expected that the number of photons is conserved down the rotational excitation ladder. The higher energy lines in the shorter wavelengths should therefore be brighter than the lower energy lines. Based on the flux of the quadruplet around ~15.3 µm, it is expected that the flux of the quadruplet at ~ 12.65 µm would be $\frac{15.3}{12.65}\approx 1.2$ times higher. This does not hold here, since these expected brighter lines are not observed. Additionally, prompt emission results in an asymmetry of the OH quadruplets, which is not found in our MIRI/MRS data (Carr & Najita 2014; Zhou et al. 2015; Zannese et al., in prep.). We therefore conclude that OH lines are not primarily excited by H₂O photodissociation. The detected OH lines are more likely excited by collision or by chemical pumping through the O + H₂ → OH + H reaction. The non-detection of OH prompt emission could indicate that small dust grains, which could be blocking the emission of the H₂O-flux from deeper layers, could also be strongly attenuating the UV radiation field.

4.6 Reduced CO₂

Similar to H₂O, the precursor for CO₂ formation is OH. However, this gas reaction is favoured over the formation of H₂O for T ≤ 250 K (Bosman et al. 2022b). In Table 1 we find the excitation temperature for CO₂ to be in this range, while the excitation temperatures for H₂O are higher in most regions. A lack of OH in the deeper and colder layers of the disk due to UV shielding of H₂O, could prevent the formation of CO₂ since it is expected to form in the deeper part of the disk.

Bosman et al. (2022b) find that the self-shielding of H₂O inhibits the formation of CO₂ due to a lack of OH, while the abundance of CO₂ can still be reduced due to dissociation deeper into the disk due to a lack of self-shielding. This is in contrast with other modelling efforts, where the CO₂ Q-branch is generally overproduced compared to the spectrum of Sz 98 (e.g. Woitke et al. 2018; Anderson et al. 2021). The temperature structure of the disk is critical: following Glassgold & Najita (2015), when Bosman et al. (2022b) add additional chemical heating in their work, the thermochemical equilibrium tips towards formation of H₂O, since this reaction is favoured for higher temperatures. These reasons combined (a change in temperature structure and self-shielding) result in a significantly reduced line flux for CO₂. As a result, the column density of the CO₂ emission is expected to be ~10¹⁶ cm⁻² (which is slightly lower than expected from other models, Anderson et al. 2021), at lower temperatures. As presented in Table 1, we indeed find these lower temperatures, although the column density is higher. However, based on the confidence in the χ² map due to the degeneracy between column density and emitting area presented in Fig. A.2, the column density could very well be lower, and perhaps in the ~10¹⁶ cm⁻² range if the emission was more extended. It is therefore not certain that this discrepancy is truly present.

While the CO₂ abundance in the inner regions of T Tauri protoplanetary disks can be greatly reduced due to the aforementioned mechanisms, Bosman et al. (2022b) still find a relatively bright Q-branch compared to the H₂O emission when only including H₂O self-shielding in their models, which does not match what we find in the disk of Sz 98. In their base model, CO₂ emits more strongly from cold gas farther out. When they assume that the emission for both species is contained within the warmer regions inside the H₂O snow line, a similar ratio of line-strengths is found as for Sz 98, indicating that the CO₂- and H₂O-emission is contained in a radially less extended region of the inner disk. Bosman et al. (2022b) attribute this to physical processes being active, for example the ‘cold-finger effect’. In their view, this would lock oxygen in the form of H₂O-ice at the midplane H₂O ice line, reducing the amount of oxygen available for CO₂ formation and limiting its extent. However, as stated previously, far-IR lines would provide better constraints related to this. Based on Table 1 for H₂O and the Q-branch-only fit in Fig. A.2 for region 2, both H₂O and CO₂ seem to be contained with a similar radial extent. However, the R value is related to an emitting area, not the location. The emission could therefore also originate from a ring further out in the disk, for example due to a small cavity. A theoretical exploration of the H₂O versus CO₂ emission variations using thermochemical disk models will be presented in Vlasblom et al. (in prep.).

Fig. 8 Close-up of the H2O and OH emission around ~13 µm. The horizontal dotted line indicates the expected flux of the OH quadruplet in case of prompt emission following H2O photodissociation by Ly-α photons.

4.7 Lack of other organics

Due to the shielding of H₂O either due to itself, or small dust, less atomic oxygen and OH are available for the formation of both CO₂ and CO. This means that the carbon is free to form hydrocarbons and other organic molecules (Duval et al. 2022). However, aside from HCN and CO₂, these species are not detected in the spectrum of Sz 98. GW Lup, as presented by Grant et al. (2023), shows detectable C₂H₂, which is lacking in the spectrum of Sz 98. After subtracting the best-fit slab models, we find an upper limit for C₂H₂ in the 13.65–13.72 µm region of ~0.3 × 10⁻¹⁴ erg s⁻¹ cm⁻² normalised to a distance of 140 pc. As can be seen in Tabone et al. (2023), this is among the lowest C₂H₂ fluxes in the sample of Banzatti et al. (2020). When detected, C₂H₂ has been fit with high temperatures in the past (Salyk et al. 2011), indicating it is located higher in the disk or at smaller radii. If the former, it is simply not abundant enough to be detected in Sz 98. On the other hand, modelling efforts expect C₂H₂ to lie deeper in the disk than HCN and CO₂ (Woitke et al. 2018). It is therefore possible that the organics are present, but blocked by the small dust deeper down in the inner disk. However, C₂H₂ was detected in the MIRI/MRS spectrum of EX Lup (Kóspál et al. 2023), despite the dust likely having a slightly smaller grain size based on the results in Kessler-Silacci et al. (2006) and Kóspál et al. (2023), indicating more opaque dust.

On the other hand, the gas in the outer disk of Sz 98 has been inferred to have a C/O > 1 (Miotello et al. 2019), whereas the inner disk may be poor in carbon to begin with. The infall of icy pebbles is expected to mainly replenish the inner disk’s oxygen abundance rather than its volatile carbon abundance (e.g. Öberg & Bergin 2021), resulting in a lower volatile elemental C/O ratio than expected from observations of the outer disk. The diversity in the C₂H₂ and HCN fluxes in Spitzer samples could be caused by differing C/O ratios between disks (Carr & Najita 2011; Walsh et al. 2015). The larger the deviation from solar C/O ratios (~0.5), the larger the range in flux ratios. The samples examined in Carr & Najita (2011) (line fluxes: H₂O, CO₂, HCN) and Salyk et al. (2011) (line fluxes: CO₂, HCN, OH, C₂H₂) have line flux ratios spanning from 0.1 to 10. However, the lack of (or at the very least extremely weak) C₂H₂ and CH₄ in Sz 98 indicates that a larger deviation from solar C/O is to be expected. According to Fig. 11 of Anderson et al. (2021), this could correspond more to a C/O of approximately 0.14 or lower, depending on the other physical properties of the inner disk. They find the HCN/H₂O flux ratio to be a good tracer for the C/O ratio in the inner disk. We use the same region of the spectrum to calculate this flux ratio, focusing on the ~13.9 µm and ~ 17.23 µm lines of HCN and H₂O, respectively. After subtracting the slab model fits presented in Fig. 2 to get a ‘clean’ HCN feature, the HCN/H₂O flux ratio is ~1.7. Assuming their fiducial model and that this is only caused by the C/O ratio, this would indeed indicate that C/O < 0.14 (Anderson et al. 2021). This is in contrast with the observations of Najita et al. (2013), where the larger disk mass should result in a larger C/O ratio. They link this to the idea that disks of higher mass could more readily form larger planetesi-mals, depleting the inner disk of gaseous H₂O. It could therefore be the case that this is not true for Sz 98.

However, the line fluxes may change when altering the fiducial model presented in Anderson et al. (2021). In order to investigate this, Anderson et al. (2021) varied a selection of parameters in their models, and we discuss some of their conclusions in relation to the observations here. For example, a larger inner gas radius of 0.5 au, similar to what has been suggested for the dust cavity of Sz 98 (van Terwisga et al. 2019), instead of 0.2 au reduces the line fluxes of most species (and would reduce the line flux of H₂O at shorter wavelengths, see e.g. Banzatti et al. 2017, 2023), but the influence on more extended species is limited (e.g. C₂H₂ and OH). A larger inner gas radius is therefore unlikely in Sz 98, since the relative line fluxes lean towards stronger H₂O and HCN instead, both of which are found to be less extended (Anderson et al. 2021). On the other hand, when assuming the gas temperature is similar to the dust temperature, Anderson et al. (2021) find a lower line flux for C₂H₂ and OH, while the other molecules are unchanged. This could also cause the larger difference in line flux ratios, allowing the C/O ratio to be slightly higher than 0.14. However, the HCN/H₂O is less sensitive to this change. On the other hand, based on ProDiMo models Antonellini et al. (2023) find the HCN/H₂O to be more sensitive to the dust opacity, where the ratio decreases for less opaque dust. While this may be part of the cause why the flux ratio is relatively low, it is still possible that the C/O is much lower than the solar value, limiting the abundances of hydrocarbons.

In Appendix D the CO snow line is estimated to lie at ~20 au, which lies in the first ring located at ~27 au separated from the inner disk by the gap at ~16 au. This gap could block CO from reaching the inner disk. It is therefore possible that H₂O-ice migrated inwards first, reaching the inner disk, sublimating, and replenishing the oxygen reservoir prior to formation of the gaps; while CO-ice got trapped in the ring before it could migrate further. This would prevent more CO from reaching the inner disk, where a larger amount of oxygen would already be present from sublimated H₂O-ice, therefore reducing the C/O ratio. The timescale over which the different ices are delivered to the inner disk compared to the timescale over which the gaps formed proves crucial to understanding the composition of the inner disk. Combining studies discussing relative abundances of gas and ice species over time such as Eistrup et al. (2018), Eistrup & Henning (2022) with gap formation could shed more light on the effects on the species available to accreting planets in the inner disk.

5 Conclusions

We presented the JWST MIRI/MRS spectrum of the inner disk of Sz 98, a disk previously only observed with the low-resolution mode of Spitzer IRS. The improved resolution and sensitivity of the MRS reveal a rich spectrum full of both ro-vibrational and pure rotational H₂O lines superposed on the continuum. Aside from H₂O, we detect CO, CO₂, HCN, and OH. Despite the disk’s large size, the thick H₂O column indicates that grains have likely drifted in towards the star, allowing H₂O ice to sublimate, forming an optically thick, potentially self-shielding layer of gaseous H₂O. Additionally, the spectrum likely probes different H₂O reservoirs from the inner parts outwards, when analysing data from shorter to longer wavelengths. This property must be considered in future work when fitting the entire MIRI/MRS spectral range. The line fluxes at longer wavelengths are still quite high, indicating that most of the H₂O is present at the surface, and could even be past the mid-plane snow line.

We find several signs pointing towards limited photodissociation of H₂O, due to self-shielding and/or dust extinction. First, its column density is larger than the few times 10¹⁷ cm⁻² required for self-shielding. Second, the H₂O line fluxes are relatively low, potentially due the presence of small dust blocking the deeper parts of the H₂O column and causing the extinction of the UV flux. Finally, the CO₂ emission is relatively weak in this disk, due to limited availability of OH for its formation, locking most of the oxygen up in H₂O.

The lack of other organic molecules, most notably C₂H₂, in the spectrum of the inner disk of Sz 98 is indicative of a low volatile elemental C/O ratio, potentially 0.14 or lower. The use of line fluxes to determine the C/O ratio has drawbacks, since they are dependent on other disk properties as well. However, a subsolar C/O ratio seems likely for the Sz 98 inner disk, while the outer disk exhibits a high C/O and low C/H. Sz 98 is not unique in this regard; many of the disks observed with Spitzer show similar H₂O and OH/H₂O line fluxes, although for the inner disk of Sz 98 the line fluxes are on the lower end compared to the sample in Banzatti et al. (2017). If organics are also lacking in the other disks in the same range, this could mean their C/O ratios are similarly low. The timescale over which ices migrate versus the timescale over which substructures form likely influences the C/O ratio of the inner disk.

Acknowledgements

D.G. thanks John H. Black for his insightful comments. 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 and #1050. The following National and International Funding Agencies funded and supported the MIRI development: NASA; ESA; Belgian Science Policy Office (BELSPO); Centre Nationale d’Études Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luftund Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y Competividad; 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. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.0.01458.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. A.C.G. has been supported by PRIN-INAF MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-generation instruments” and from PRIN-INAF 2019 “Spectroscopically tracing the disk dispersal evolution (STRADE)”. G.B. thanks the Deutsche Forschungsgemeinschaft (DFG) – grant 138 325594231, FOR 2634/2. E.v.D. acknowledges support from the ERC grant 101019751 MOLDISK and the Danish National Research Foundation through the Center of Excellence “InterCat” (DNRF150). D.G. would like to thank the Research Foundation Flanders for co-financing the present research (grant number V435622N). 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. I.K., A.M.A., and E.v.D. acknowledge support from grant TOP-1 614.001.751 from the Dutch Research Council (NWO). I.K. and J.K. acknowledge funding from H2020-MSCA-ITN-2019, grant no. 860470 (CHAMELEON). B.T. is a Laureate of the Paris Région fellowship program, which is supported by the Île-de-France Region and has received funding under the Horizon 2020 innovation framework program and Marie Sklodowska-Curie grant agreement No. 945298. O.A. and V.C. acknowledge funding from the Belgian F.R.S.-FNRS. I.A., D.G. and B.V. thank 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). L.C. acknowledges support by grant PIB2021-127718NB-I00, from the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033. T.P.R acknowledges support from ERC grant 743029 EASY. D.R.L. acknowledges support from Science Foundation Ireland (grant number 21/PATH-S/9339). D.B. has been funded by Spanish MCIN/AEI/10.13039/501100011033 grants PID2019-107061GB-C61 and No. MDM-2017-0737. M.T. acknowledges support from the ERC grant 101019751 MOLDISK

Appendix A χ² maps of fits

The χ² maps per molecule and per region are presented here. H₂O across all regions is presented in Fig. A.1. The confidence of HCN and CO₂ detected in region 2 can be found in Fig. A.2.

Fig. A.1 χ2 plots of H2O over the different regions of the disk (top left: region 1, top right: region 2, bottom left: region 3, bottom right: region 4). The red, yellow, and green lines indicate the 1σ, 2σ, and 3σ confidence contours, respectively. The white contours show the emitting radii in astronomical units (0.1 to 2.0 au). The black cross corresponds to the best fit. The black dotted line corresponds to the upper limit of the column density based on the H218O and H216O flux ratio as described in Sect. 3.1.

Fig. A.2 χ2 plots of HCN and CO2 detected in region 2. The red, yellow, and green lines indicate the 1σ, 2σ, and 3σ confidence contours, respectively. The white contours show the emitting radii in astronomical units (0.1 to 2.0 au). The black cross corresponds to the best fit. The bottom plot shows the map for CO2 when fitting the Q-branch only.

Appendix B Masked features

The wavelengths at which spurious data reduction artefacts have been masked are: [5.0091,5.01071]; [5.018,5.019]; [5.112,5.15]; [5.2157,5.2184]; [5.2267,5.2290]; [5.2441,5.2471]; [5.2947,5.2974]; [5.3742,5.3777]; [5.3836,5.3877]; [5.4181,5.4210]; [5.5644,5.5674]; [5.5925,5.5966]; [5.8252,5.8267]; [5.8669,5.8689]; [5.9,5.916]; [5.9282,5.9314]; [5.9691,5.9728]; [6.0357,6.0394]; [6.0430,6.0462]; [6.1012,6.1044]; [6.1311,6.1421]; [6.26,6.31]; [6.3740,6.3757]; [6.3783,6.3810]; [18.8055,18.8145]; [19.004,19.012]; [21.974,21.985]; [25.69824,25.71313].

Due to the calibration files being based on an A star, some of its features are propagated into the target here. The majority of these features are present in the short wavelengths, resulting in most of the blanked lines being present here.

Appendix C CO₂ column density

The column density of CO₂ is difficult to constrain due to the faintness of the features and the overlap with H₂O. In Fig. C.1 it can be seen that distinguishing between the optically thin and thick cases is non-trivial.

Fig. C.1 Comparison between the optically thick and thin cases for CO2 and 13CO2 assuming a ratio of 68, overlain on the data with H2O, HCN, and OH subtracted. Due to the faintness of the features and blending with H2O lines it becomes difficult to constrain the best model.

Appendix D Sz 98 from ALMA data

The continuum image shown in Fig. 6 has been created using the continuum spectral window (centred at 232.984 GHz, 1286.75 µm) contained in the ALMA archival data set 2018.1.01458.S (PI: Hsi-Wei Jen). The data set was calibrated using the provided pipeline scripts and the specified Common Astronomy Software Applications (CASA) (McMullin et al. 2007). The continuum image was created using CASA version 6.5.2.26 and the TCLEAN-task, where we have used the Briggs weighting scheme and a robust parameter of +1.0. We cleaned down to a threshold of ~ 1 × the RMS in the initial, dirty image. The inclination of the disk is 47.1° with a position angle of 111.6° (Tazzari et al. 2017).

The cleaned image has a resolving beam of 0.079″ × 0.063″ (-27.135°). The emission peaks at 2.167 mJy beam⁻¹ and the image has a root mean square of 0.071 mJy beam⁻¹. We have estimated the locations of the rings and gaps (displayed in Fig. 6) using a deprojected, azimuthally averaged radial profile. The radial profile was created using bin sizes of half the width of the beam’s minor axis. As displayed in Fig. 6, we infer the millimetre emission to consist of an inner compact core, surrounded by two rings at approximately 27 and 103 au, respectively. The gaps are expected to be located at approximately 16 and 88 au.

The location of the H₂O, CO₂, and CO snow lines are estimated from a radiative transfer model that fits the spectral energy distribution, mm dust radial profiles and the JWST continuum. Assuming that these molecules freeze-out at temperatures of, respectively, 150 K, 72 K, and 20 K (e.g. Collings et al. 2004), the model suggests that the H₂O snowline is located between approximately 1.05 and 1.11 au, the CO₂ snowline between approximately 2.17 and 2.29 au, and the CO snowline between approximately 20.64 and 21.40 au.

References

All Tables

All Figures

	Fig. 1 Full MIRI/MRS spectrum of the Sz 98 disk. The extent per region is indicated, along with the detected species. The Spitzer LRS spectrum from the CASSIS database is shown in light grey (Lebouteiller et al. 2011).
In the text

Fig. 2 Slab fits for the four different wavelength ranges. The top and bottom panels per region are the data and total model, and the individual slab models; respectively. Spurious spikes and features from the data reduction have been blanked out (see the list in Appendix B). The horizontal lines indicate the regions used to fit the molecules, or the region where σ is estimated after subtracting the slab models.

In the text

	Fig. 3 Best-fit temperature and emitting radius for H2O per region of the inner disk. The error bars indicate the 1σ contours of the χ2 maps.
In the text

Fig. 4 Flux ratios of H216O lines of the same upper level in the data compared to slab models (top); and H216O/H218O ratio in the data and slab models (bottom). The properties of the transitions can be found in Table 2. The dotted black vertical line indicates the upper limit of the column density based on the 1σ level of the spectrum and the nondetection of the H218O line. The x-axis shows the column density of H216O. In the bottom panel, H216O/H218O ratio of 550 is assumed.

In the text

Fig. 5 Different regions where H218O is expected to show lines sufficiently isolated from H2O. These lines are not seen in the spectrum, compared a H218O slab models assuming the same best-fit parameters as H2O in these regions, assuming a 16O/18O of 550. The models have been multiplied by 5 for visibility. The H2O slab models are based on the best-fit parameters presented in Table 1 per wavelength region.

In the text

	Fig. 6 ALMA millimetre dust continuum emission from Sz 98 centred at 232.984 GHz (1286.75 µm). The dotted lines indicate gaps (G), the solid lines indicate rings (R), and their radii are included in au. Appendix D describes how the image is created.
In the text

	Fig. 7 Line fluxes of H2O (top) and line flux ratio OH/H2O (bottom) and accretion luminosities from the disk sample in Banzatti et al. (2017) compared to Sz 98 (red dot). The size of the circle is proportional to the stellar mass, and the arrows indicate upper and lower limits. All fluxes have been normalised to 140 pc.
In the text

	Fig. 8 Close-up of the H2O and OH emission around ~13 µm. The horizontal dotted line indicates the expected flux of the OH quadruplet in case of prompt emission following H2O photodissociation by Ly-α photons.
In the text

Fig. A.1 χ2 plots of H2O over the different regions of the disk (top left: region 1, top right: region 2, bottom left: region 3, bottom right: region 4). The red, yellow, and green lines indicate the 1σ, 2σ, and 3σ confidence contours, respectively. The white contours show the emitting radii in astronomical units (0.1 to 2.0 au). The black cross corresponds to the best fit. The black dotted line corresponds to the upper limit of the column density based on the H218O and H216O flux ratio as described in Sect. 3.1.

In the text

	Fig. A.2 χ2 plots of HCN and CO2 detected in region 2. The red, yellow, and green lines indicate the 1σ, 2σ, and 3σ confidence contours, respectively. The white contours show the emitting radii in astronomical units (0.1 to 2.0 au). The black cross corresponds to the best fit. The bottom plot shows the map for CO2 when fitting the Q-branch only.
In the text

	Fig. C.1 Comparison between the optically thick and thin cases for CO2 and 13CO2 assuming a ratio of 68, overlain on the data with H2O, HCN, and OH subtracted. Due to the faintness of the features and blending with H2O lines it becomes difficult to constrain the best model.
In the text