JWST-MIRI Spectroscopy of Warm Molecular Emission and Variability in the AS 209 Disk

We present MIRI Medium-resolution Spectrograph observations of the large, multi-gapped protoplanetary disk around the T Tauri star AS 209. The observations reveal hundreds of water vapor lines from 4.9–25.5 μm toward the inner ∼1 au in the disk, including the first detection of rovibrational water emission in this disk. The spectrum is dominated by hot (∼800 K) water vapor and OH gas, with only marginal detections of CO2, HCN, and a possible colder water vapor component. Using slab models with a detailed treatment of opacities and line overlap, we retrieve the column density, emitting area, and excitation temperature of water vapor and OH, and provide upper limits for the observable mass of other molecules. Compared to MIRI spectra of other T Tauri disks, the inner disk of AS 209 does not appear to be atypically depleted in CO2 nor HCN. Based on Spitzer Infrared Spectrograph observations, we further find evidence for molecular emission variability over a 10 yr baseline. Water, OH, and CO2 line luminosities have decreased by factors of 2–4 in the new MIRI epoch, yet there are minimal continuum emission variations. The origin of this variability is yet to be understood.


INTRODUCTION
The presence of water in circumstellar disks plays a fundamental role in the outcome of planet formation.Water is expected to be a major oxygen carrier in disks-therefore key to set the local gas-phase and solid C/O and N/O ratios in different disk regions (e.g.Herbst & van Dishoeck 2009;van Dishoeck et al. 2014;Öberg & Bergin 2021)-and undeniably an essential ingredient for setting planetary habitability.In the outer disk, beyond the water snowline, the stickiness of water ice (among other factors) is thought to facilitate the growth of pebbles and planetesimals (Dominik & Tielens 1997;Birnstiel et al. 2010;Gundlach & Blum 2015).Moreover, the migration efficiency of these bodies might be regulated by the local concentration of water ice (Stevenson & Lunine 1988;Ciesla & Cuzzi 2006).This solid water reservoir, however, is not observable in most disks.Instead, water can be most readily characterized via its gas-phase emission spectrum, which arises from the warm, innermost disk regions (within 10 au) in the mid-IR (Pontoppidan & Blevins 2014), and from the outer disk at far-IR wavelengths (Hogerheijde et al. 2011;Du et al. 2017).
Our current view of inner-disk water comes primarily from the Infrared Spectrograph onboard the Spitzer Space Telescope (Spitzer IRS ) (Houck et al. 2004), which revealed that bright, mid-IR rotational water emission is ubiquitous among low-mass T Tauri stars (Carr & Najita 2008;Pontoppidan et al. 2010b).In general, the intensity of mid-IR water emission lines in these sources is best reproduced by column densities of order 10 18 cm −2 , warm excitation temperatures ranging from T ex ≈ 200 to 800 K, and emitting areas in the disk of around 10 au 2 , corresponding to outer circular radii of ≲ 2 au (Carr & Najita 2011;Salyk et al. 2011;Najita et al. 2011;Antonellini et al. 2015).These emitting areas tend to be consistent with line kinematic information inferred from high-resolution ground-based spectra (Pontoppidan et al. 2010a;Banzatti et al. 2014;Najita et al. 2018;Salyk et al. 2019;Banzatti et al. 2023a).Mid-IR observations have further revealed that the intensity of water emission is anti-correlated to the mass and evolutionary stage of the central star (Pontoppidan et al. 2010b;Salyk et al. 2011;Fedele et al. 2011).Determining the abundance and distribution of water vapor * NASA Hubble Fellowship Program Sagan Fellow in disks can thus provide substantial insights regarding the disk's dynamical and thermal evolution.
The characterization of water vapor and the trends described above have so far encountered two critical constraints.First, due to the low spectral resolution and sensitivity of Spitzer IRS, most of the water lines observed with this instrument are blended, and model fits were dominated by the brightest, most optically thick water lines that emit near the disk surface (e.g.Salyk et al. 2011;Carr & Najita 2011).The end result is that slab models infer 10 − 100× lower abundances relative to the total water column predicted by thermochemical simulations that account for chemical heating and H 2 O self-shielding (Meijerink et al. 2009;Blevins et al. 2016;Woitke et al. 2019;Bosman et al. 2022).
Second, a complete understanding of water vapor in disks has been precluded by the limited wavelength coverage of Spitzer IRS (Banzatti et al. 2023a).While Spitzer IRS could readily detect bright rotational lines tracing lukewarm (400-600 K) water vapor in disks, a hotter (T ex ∼ 1000 K), more compact water reservoir is only revealed by ro-vibrational water emission within 2-8 µm (Carr et al. 2004;Salyk et al. 2008;Mandell et al. 2012), which had been only accessible to groundbased spectrographs.Similarly, a cold water reservoir (T ex ≈ 200 K) associated with sublimation near the snowline and photochemistry in the outer regions of some disks was mainly observable with far-IR facilities like Herschel (Zhang et al. 2013;Du et al. 2017).No single instrument until now has enabled a comprehensive analysis of distinct water emission reservoirs with sufficient sensitivity and spectral resolution.
Thanks to its unparalleled capabilities, the James Webb Space Telescope (JWST) has the potential to overcome the difficulties mentioned above and accurately characterize excitation conditions of water vapor in disks.In particular, the Mid-Infrared Instrument (MIRI; Rieke et al. 2015) Medium-resolution Spectrometer (MRS; Wells et al. 2015) offers medium spectral resolution (R = 1500 − 3700) from 4.9 to 28 µm, allowing the simultaneous analysis of ro-vibrational (bending mode) and rotational water lines spanning a vast range of upper energy levels (E u = 10 3 − 10 4 K) for the first time.Similarly, the sensitivity of MIRI provides the opportunity to identify and characterize multiple optically thin water lines, and hence trace water closer to the planet-forming disk midplane.Indeed, JWST MIRI has already started to provide robust evidence for a lowenergy water line excess in compact disks (Banzatti et al. 2023b) as well as the first detections of weak isotopologue lines (Grant et al. 2022), small organics (Kóspál et al. 2023), and ionization tracers (Berné et al. 2023) in the terrestrial-planet forming region of disks.
In this paper, we model the emission of H 2 O, OH radical (hereafter OH), CO 2 , C 2 H 2 , and HCN-with particular focus on water vapor and OH-using MIRI MRS data of the inner disk of AS 209.The disk of AS 209 is the first of four targets observed as part of the JWST follow-up to the ALMA MAPS Large Program ( Öberg et al. 2021), and the only one with prior inner-disk chemistry constraints from Spitzer IRS in this sample.In Section 2, we describe the observations and data reduction process.The spectral modeling technique and retrieval framework are detailed in Section 3, and Section 4 presents the excitation analysis results.In Section 5 we explore whether the mid-IR spectrum of AS 209 differs from other T-Tauri disks, and compare our findings with prior Spitzer IRS observations.We find evidence for significant variability in molecular emission features, and present some possible explanations for this anomaly.

The Disk of AS 209
The low-mass (0.96M ⊙ ; Fang et al. 2018) T Tauri star AS 209, located at 121 pc (Gaia Collaboration et al. 2018), hosts an unusually bright and extended protoplanetary disk, with a mm-dust radius of R dust = 139 au and a CO 2 − 1 gas disk that extends out to 272 au (Andrews et al. 2018;Huang et al. 2018;Guzmán et al. 2018).Based on ALMA observations, the disk exhibits a least seven deep dust gaps from 9 out to 137 au (Huang et al. 2018;Sierra et al. 2021).Many molecules have been detected in the AS 209 outer disk (beyond 10 au) at sub-millimiter wavelengths (e.g.Öberg et al. 2011;Huang et al. 2017;Bergner et al. 2018;Öberg & Bergin 2021), with most exhibiting multiple ring-like emission substructures when observed at high spatial resolution (Law et al. 2021).
There is evidence that at least one of the dust gaps is associated with ongoing planet formation at 100 au (Favre et al. 2019), and a candidate circumplanetary disk has been detected in molecular gas at ∼ 200 au (Bae et al. 2022).The possibility that additional planets may be forming within the innermost ∼ 10 au makes AS 209 an attractive target for IR follow-up observations.Previous Spitzer IRS observations of this disk suggested that water vapor is abundant in the innermost 2 au of the disk, while OH, CO 2 , and HCN gas were only tentatively identified (Salyk et al. 2011).Far-IR observations with Herschel, on the other hand, have only placed an upper limit of 75.3 mK km/s for the in-tegrated flux of the 1 1,0 − 1 0,1 water vapor line in the far-IR, and ro-vibrational lines were not detected either (Banzatti et al. 2023a).
The data was reduced using version 1.11.0 (CRDS context jwst 1105.pmap) of the JWST pipeline (Bushouse et al. 2023).For further details regarding the data reduction and 1D spectra extraction, we refer the reader to Pontoppidan et al. (2023).Being an atypically bright point-source as observed by MIRI MRS, the spectrum of AS 209 is severely affected by residual fringes, even after applying the fringe correction procedure ResidualFringeStep built into the JWST pipeline.We observe broad, quasi-periodic artifacts in all channels, and significant high-frequency variations in brightness versus wavelength around 10 − 12 µm and in all of channel 4.These artifacts complicate the detection of real emission lines considerably, as well as the location of the true continuum baseline in each channel.We thus apply an empirical fringe correction using two asteroid calibrators from JWST GO program 1549 (P.I.Klaus Pontoppidan), which is also fully described by Pontoppidan et al. (2023).
The top panel of Figure 1 shows the full MIRI MRS spectrum of the AS 209 disk, and the bottom panel shows the continuum-subtracted spectrum.To subtract the continuum, we developed an iterative Gaussian Processes decomposition technique that allows us to fit a baseline to the data without having to identify all linefree regions amid the noisy background.The details are presented in Appendix A. Even after the correction for residual fringes, some regions of the spectrum, including most of channel 4, still exhibit multiple artifacts, bad pixels, and overall low SNR.We decide to exclude the 8. 5 -12.3, 17.5 -20.0, and 25.5 -28.1 µm wavelength ranges (shaded regions in Figure 1) from our analysis.
Figure 2 provides a closer view at the regions of the MIRI MRS spectrum analyzed in this paper.In each panel, we include a template water vapor model (in violet) with an excitation temperature of 700 K and column density of 10 18 cm −2 , scaled 'by-eye' to match the data.These values were chosen based on fits to highenergy water lines in other disks observed with MIRI MRS (Banzatti et al. 2023b).As shown in the top panel, MIRI MRS reveals the presence of ro-vibrational water emission in the disk of AS 209 for the first time, which is too weak to be detected from the ground (Banzatti et al. 2023a).We label some of the most prominent emission lines from other species that could be visually identified, namely CO, H 2 , OH, and a few atomic lines.Based on visual inspection alone, we find no evidence of significant CO 2 , HCN, nor hydrocarbon emission lines in the spectrum, in line with the observations provided by Spitzer IRS.Since the 4.9 -8.5 µm region is dominated by the ro-vibrational bending mode of water, hereafter we refer to this wavelength range as the ro-vibrational region.The 12.3 -17.5 and 20.0 -25.5 µm wavelength ranges are dominated by rotational water emission, and we refer to these as the rotational region.

ANALYSIS
In the following we analyze the spectrum in detail to derive the column density (N mol ), excitation temperature (T ex ), and emitting area in the disk (A d ) of the two species found based on visual inspection: H 2 O and OH.We then perform a deeper search for emission of other molecules of interest: H 2 , CO 2 , C 2 H 2 , and HCN, and estimate their column densities when possible.

Line Modeling
To simulate the molecular emission, we developed and deployed the GPU-accelerated LTE slab modeling Python package iris (Munoz-Romero et al. 2023).The line modeling prescription follows that presented by Tabone et al. (2023), adopting a detailed treatment of optical depth effects and line overlap and using spectroscopic information from the HITRAN database (Gordon et al. 2022).In brief, we generate a high-resolution (R ∼ 10 5 ) opacity-weighted flux model for each species, with an optical-depth grid described by for each line n, where λ 0,n is the line center wavelength and ∆V is the line full width at half maximum (FWHM).Here, τ 0,n is the line center optical depth with where A u,l is the transition Einstein coefficient, x u and x l are the upper-and lower-level populations, respectively, and g u and g l their degeneracies.The line FWHM is assumed to be the sum in quadrature of a fixed turbulent component (v turb ) and a thermal broadening component where µ is the molecular weight of each species, k B is the Boltzmann constant, and m p the proton mass.As the observations are spectrally unresolved, the turbulent width cannot be directly determined, and we decide to use the same value assumed by Salyk et al. (2011) and Tabone et al. (2023), v turb = 4.71 km/s (2 km/s standard deviation of the Gaussian line profile).Since most of the lines are expected to be optically thick, we caution that the derived column densities scale roughly as 1/∆V .Finally, the flux model is calculated as where B ν (T ex ) is the full Planck function and d the distance to AS 209, assumed to be 121 pc.The above prescription is computationally expensive-each species adds hundreds to thousands of transitions, all of which are modeled at high-resolution as described above-yet necessary given the plethora of emission lines that overlap within their intrinsic line widths.Ignoring these effects results in a severe overestimation of peak line intensities, even at moderately high column densities (≳ 10 18 cm −2 ).We convolve the emission models with Gaussian kernels to approximate the resolving power of MIRI MRS (R = 3400, 2300, 1700 for the 4.9 -8.5, 12.3 -17.5, and 20 -25.5 µm wavelength ranges, respectively) and downsample the resulting models to the MIRI MRS wavelength grid (while conserving the total flux) using the Python package SpectRes (Carnall 2021) to compare with the data.

Bayesian Retrieval
We model the ro-vibrational and rotational regions independently and include different species in each wavelength range of interest.For the rotational region, we include H 2 O, OH, CO 2 , C 2 H 2 , and HCN in the 12.3 -  ).We emphasize that when multiple species (or multiple components for a single species) are included, all of them are modeled simultaneously, as modeling and subtracting each one individually would not correctly account for line opacity overlap along the line of sight.The rest of this section elaborates on the specific fitting details.
We perform Bayesian inference in a hierarchical framework, using the Python Nested Sampling package dynesty (Speagle 2020) to explore the posterior distribution of each parameter of interest.By iteratively adding molecules to the fit, we are able to better guide the sampler and avoid adding complexity to the model without sufficient evidence.To summarize the method described below, we add molecules to each emission model in succession, and each time a new molecule is included, we use the posterior distribution of the previous (simpler) fit to inform the prior distributions used in the new model.

Water Vapor
For the rotational region, we begin by fitting an H 2 O model using the following Uniform prior distributions: log 10 T ex = U(2.17,3.30) K log 10 N mol = U(13, 22) cm −2 with bounds informed by the previous observations and modeling efforts of water vapor emission in disks as discussed in Section 1.The same priors are used to model the H 2 O emission in the rotational region.Note that the temperature bounds above correspond to 150 and 2000 K, and the log area bounds correspond to equivalent radii of 0.017 -17.8 au assuming circular emitting areas.The nested sampling is performed using multiple ellipsoid decomposition and random walk sampling with n live = 10 × n dim chains (i.e.live points), where n dim is the number of parameters in the model.We use a Normal likelihood assuming a constant 20 mJy uncertainty for all pixels1 .The run is terminated once the estimated change in the remaining evidence ∆ log Z satisfies ∆ log Z ≤ 0.001(n live − 1) + 0.01.

OH
The next molecule added to the rotational region model is OH, which is also visually identified in the spectrum.To perform this fit, we begin by using the same Uniform prior distributions described by Equation 6for OH, and Normal prior distributions for H 2 O.In particular, we take the median and 3× the standard deviation of the posterior distribution of the previous H 2 O fit to construct the new Normal priors.We choose to use 3× the standard deviation of the previous results to avoid under-estimating the parameter uncertainties in the new fits.The number of live points is increased and the run is terminated using the same criteria as described above.
We find that the above fit produces multi-modal posterior distributions for OH, with three main solutions: one with T ex of ∼ 400 K, one near 700 K, and a third one at approximately 1500 K.The first two solutions have similar emitting areas of ∼ 1 au 2 and column densities of order 10 18 cm −2 , while the third one requires a column density < 10 16 cm −2 .This is not unexpected, since the gas responsible for inner disk OH emission is thought to be in non-LTE, and slab models with multiple excitation temperature components are typically required to reproduce the observed line ratios (e.g.Najita et al. 2010;Banzatti et al. 2012;Carr & Najita 2014).Of the three solutions, we focus on the one with T ex ≈ 700 K (see Section 4) for two main reasons: First, this solution has the highest likelihood and can be considered to be the global minimum in the parameter space allowed by the priors.Second, we find that the low temperature solution misses the high-energy OH lines seen at 16.8 µm, while the high temperature solution cannot reproduce the lower energy lines beyond 20 µm.The 700 K solution on the other hand provides the best overall fit to the data.We also experimented with fitting multiple OH components simultaneously, but the models did not converge.
Based on the findings described above we repeat the H 2 O + OH fit, this time using a Normal T ex prior for OH centered at 700 K with standard deviation of 200 K, and the same N mol and A d Uniform priors.With the more restrictive temperature prior we no longer recover the multi-modal solutions for OH, and thus the posterior distribution provides a better estimate of the uncertainty for the main solution.

A Cold Water Vapor Component
The fit described before results in a hot and compact water emission component, as further detailed in Section 4.1.After careful inspection of the residuals, we find minor excess emission features that may be associated with low-energy H 2 O lines around 21.2, 21.8, 23.8, and 25.5 µm, which are not captured by the best fit water model.As mentioned in Section 1, mid-IR rotational water vapor spectra in disks tends to exhibit excitation temperature gradients, sometimes associated with multiple reservoirs produced by distinct water vapor formation mechanisms.We thus explore whether there could be a colder water reservoir contributing to the spectrum.
We fix the temperature of this second component to 400 K, as is observed in other T-Tauri systems (Banzatti et al. 2023b), and fit only for the emitting area and column density using the same uniform priors as described above.We also experiment by allowing the excitation temperature to vary, and get a solution with T ex ∼ 200 K.However, given the low quality of the data in channel 4, we decide not to derive any conclusions regarding the temperature of cold water vapor from the long wavelength emission, and only use the fit to obtain an upper limit for the amount of cold water vapor.Similarly, we caution against interpreting the best fit cold H 2 O emit-ting area as a snowline location tracer, given that it is strongly correlated to the choice of T ex .

Other Molecules
For the rest of the species included in the 12.3 -17.5 µm range, CO 2 , C 2 H 2 , and HCN, we are unable to fit all N mol , T ex , A d , likely due to a mix of intrinsically low fluxes, few observable lines, and poor SNR.Specifically, their column densities and emitting areas are degenerate, and we cannot place any meaningful constraints on their excitation temperatures.Thus, we attempt to estimate only the column density of each species, using the following uniform prior: We use the same area and excitation temperature obtained for the hot water vapor component to estimate N mol for CO 2 , HCN, and C 2 H 2 , which are added to the model in that order.Note that when adding each of these species we no longer let the parameters of the previous models vary, but rather fix them using the median of their final posterior distribution.At the expected abundances, these species emit most prominently in non-overlapping wavelength regions (see Figure 4), and thus should not affect the retrieved properties of each other, nor those of H 2 O and OH.

Rotational Region
As described in Section 3, we attempt to fit the emission of hot and cold H 2 O, OH, CO 2 , HCN, and C 2 H 2 .Since we can only visually confirm the presence of hot water vapor and OH in the spectrum, we emphasize that CO 2 , C 2 H 2 , HCN, and cold H 2 O can be considered, at best, marginally detected.As a way to decide whether a species is detected, marginally detected, or non-detected, we use the Bayesian Information Criterion (BIC), which penalizes the increased complexity of each model even if it has a higher likelihood: where n dim is the dimensionality of the problem, k is the total number of observations, and L ⋆ the log-likelihood of the model evaluated using the median of the posterior distribution for each parameter.Specifically, hot water vapor is considered an unambiguous detection, and the rest of the molecules are: 1. Detected : if we can identify multiple lines, the strongest component lies at the > 3σ level after subtracting hot H 2 O from the data, and the BIC decreases by including the species in the model.2. Marginally detected : if it cannot be unambiguously identified in the data, but the BIC decreases by including it in the model.
3. Non-detected : if it cannot be unambiguously identified in the data, and the BIC does not decrease by including it in the model.
Figure 3 shows how the BIC is modified by adding additional species to the model, starting from hot water vapor.OH is detected, CO 2 , HCN, and cold H 2 O are marginally detected, and C 2 H 2 is not detected.CO and H 2 are considered detections too, but are not included in the excitation analysis.
Table 1 lists the best-fitting values and confidence intervals for each parameter of interest in the final fit.The best-fit values are chosen to be the median of the posterior distributions, while the uncertainties correspond to the 0.0015 and 0.9985 percentiles (3σ confidence intervals).We caution that the reported uncertainties reflect the assumed flux uncertainty of 20 mJy, but likely under-estimate the true parameter uncertainties.This is because the likelihood function and retrieval framework used in this work can only include Normally-distributed noise, but does not account for correlated noise nor the uncertainty related to both the absolute flux calibration and continuum baseline.Based on visual inspection of the posteriors, we do not notice any strong correlations between parameters.From here onward, the emission models referred to as 'best-fits' correspond to those evaluated using the median of the posterior distribution for each parameter.Since CO 2 and organics are only marginally/not detected, we consider their best-fit column density to be upper limits.In the case of cold H 2 O for which we can fit both emitting area and column density, we consider the observable mass (product of best-fit column density and emitting area) to be an upper limit as well.Table 1 further lists the inferred total observable mass of each molecule, and their ratio relative to hot H 2 O.
Figure 4 shows the best-fit model in the 12.3 -17.5 µm region.For clarity, the spectrum and model have been broken up into four vertically-stacked panels.Besides the best fit model, we show the individual contribution from each species as shaded areas of different colors, and provide additional panels with a zoomed-in view of the main organic bands.We find that the best fit hot H 2 O model does an excellent job at reproducing most of the emission observed in this spectral region, except in the noisier 12.3 -13.6 µm region.Similarly, the 700 K OH model can reproduce the strongest high-energy OH lines near 16.8 µm and even those blended with water near 16.0 µm.CO 2 and the organics, if present, are completely blended with water.Yet we note that adding HCN and CO 2 , which are marginally detected, improve the quality of the overall fit.
Finally, Figure 5 shows the best-fit model in the 20.0 -25.5 µm region.In this region, notice that hot H 2 O has only minor contributions, and the data is dominated by noise.We can yet observe that the OH model matches the brightest lines near 21.1 and 21.5 µm well, but under-predicts the observed fluxes of lower-energy lines beyond 23 µm by about 50%.Even with the higher noise, we notice that several bright features can be partially reproduced by the marginally-detected cold H 2 O component, particularly at 21.4, 21.85, between 23.8 -24.0, and between 25.0 -25.4 µm.

Ro-vibrational Region
Here we present the results from fitting the water vapor emission in the ro-vibrational region (4.9 -8.5 µm range).However, we emphasize that emission from the ro-vibrational bending mode of water is generally found to not be in local thermodynamic equilibrium (LTE)-hence why it is fit separately from the rest of the spectrum-and the derived properties using LTE slabs may not be good estimates of the true gas temperature, column density, and extent.In fact, it is likely that the ro-vibrational water emission traces a similar, perhaps slightly hotter water vapor reservoir seen in the rotational region-but is simply suppressed due to excitation effects (Banzatti et al. 2023a;Bosman et al. 2022).
Figure 6 shows the best-fit model in the ro-vibrational region, both compared to the data and individually.We find that the ro-vibrational water vapor emission component is best modelled by a slab with area of 9 × 10 −3   3.9 +1.1 −0.8 × 10 18 9.1 +1.5 −0.9 × 10 −3 --  au 2 , a higher excitation temperature of 1230 K, and a lower column density compared to that traced by the rotational component, of 4 × 10 18 cm −2 .Note that while this H 2 O model appears to be a good fit to the data, it severely under-predicts the observed water line fluxes in the rotational region.Similarly, the hot and cold water components fit to the rotational regions over-predict the observed water spectrum between 4.9 and 8.5 µm.The two regions cannot be reproduced by a single slab or sum of slab models.Based on the best-fit column densities of the other species in the rotational region, we do not expect any of them to be detected between 4.9 and 8.5 µm.

The Mid-IR Spectrum of AS 209 in Context
To guide the interpretation of our results, we carry out a detailed comparison between the MIRI MRS spectrum of AS 209 and those of two other systems: CI Tau and GK Tau.These provide some of the highest quality MIRI MRS data published for T-Tauri systems to date, and let us explore 1) if the properties of the hot water vapor reservoir in AS 209 are atypical, and 2) whether the weaker molecular features tentatively found in the AS 209 spectrum are consistent, within the noise, with those detected in other disk spectra.Furthermore, we choose these two systems to contrast the H 2 O and OH excitation conditions in AS 209 to those of a similarly extended and structured disk (CI Tau, with sub-mm dust radius R dust = 190 au) and a significantly more compact disk (GK Tau, R dust = 13 au) (Long et al. 2019).
The GK Tau and CI Tau MIRI MRS spectra were obtained in program JWST Cycle 1 program GO-1640 (P.I.Andrea Banzatti), and were recently presented by Banzatti et al. (2023b) in the context of water delivery.The spectra were reduced and de-fringed using the same calibration pipeline applied in this work to AS 209 (Pontoppidan et al. 2023).To accurately compare all three spectra, we first scale the GK Tau and CI Tau data to the same distance as AS 209 (121 pc), assuming a distance of 160 pc for CI Tau and 129 pc for GK Tau (Gaia Collaboration et al. 2018).Next, we use the bright and optically-thick high-energy H 2 O lines near 16 µm to remove the effects of differing line emitting areas (Banzatti et al. 2020).Specifically, we take the 16.05 -16.6 and 17.0 -17.3 µm spectral ranges, and find the scaling factor that minimizes the sum of squared differences between each spectrum and that of AS 209.CI Tau is scaled by a factor of 0.44, while the GK Tau spectrum does not require any additional scaling.
A first test is to determine whether the scaling factors necessary to match the optically-thick H 2 O line luminosities are consistent with the expectation based on the stellar accretion luminosity (L acc ) of each system.In particular, Banzatti et al. (2023b) find that the optically-thick H 2 O log line fluxes follow a relatively tight linear correlation with log 10 L acc , with slope ∼ 0.59. Figure 7 shows the observed and expected scaling factors, assuming log 10 L acc values of -1.12, -0.7, and -1.38 for AS 209, CI Tau, and GK Tau, respectively, and a typical 20% uncertainty in log 10 L acc (Fang et al. 2018).While the observed scaling factors are slightly lower than expected from L acc effects alone, they are still fully consistent within the uncertainty.That is, we find no evidence that the hot water vapor reservoir of AS 209 is significantly more compact or diffuse relative to these CTTs.
Note however that the retrieved emitting area and column density of hot H 2 O in AS 209 are roughly an order of magnitude smaller and higher, respectively, compared to those reported for GK Tau and CI Tau by Banzatti et al. (2023b).We have reanalyzed the spectra of the latter two using the updated line modeling framework presented in this paper, and conclude that the difference does not reflect a real discrepancy in the water vapor distribution, but rather arises from a) a different intrinsic line FWHM used (1 km/s v. 4.7 km/s used in this work), and b) the distinct treatment of line overlap and opacities.These new water vapor retrievals will be presented in a forthcoming paper (Muñoz-Romero et al. in prep.).
Figure 8 shows the comparison between the MIRI MRS spectrum of AS 209 and the scaled spectra of CI Tau and GK Tau, from 13.6 to 17 µm.For both disks, we find that their spectra are remarkably similar to that of AS 209.The OH, HCN and CO 2 features detected in both CI Tau and GK Tau are consistent with the AS 209 data.In fact, based on visual inspection alone, the main HCN and CO 2 branches seem to be slightly brighter in AS 209 compared to GK Tau.We surmise that neither species is atypically depleted in AS 209, but rather they are only marginally detected due to the lower data SNR.On the other hand, it appears that C 2 H 2 may indeed be suppressed in AS 209, as indicated by our best fit model as well, but deeper observations will be necessary to confirm this and any possible implications regarding the gas phase C/O ratio in the inner disk.

Emission Variability
We conclude by comparing the MIRI MRS spectrum of AS 209 and best fit models to the previous observations of this disk with Spitzer IRS.AS 209 was observed as part of the cores to disks (c2d) Spitzer Legacy program (Evans et al. 2003), providing coverage from 9 to 37 µm.The inner disk molecular emission was subsequently presented and characterized by Pontoppidan et al. (2010b) and Salyk et al. (2011).In particular, Salyk et al. (2011) reported the detection of water vapor and CO gas in AS 209, as well as a tentative detection of CO 2 .For H 2 O, Salyk et al. (2011) reported a best-fit column density of 10 20.9 cm −2 , an excitation temperature of 250 K, and an emitting area of 9 au 2 , which was obtained by fitting the integrated fluxes of 65 individual water lines.We find that this model strongly over-predicts the water line emission seen in the MIRI MRS spectrum.Previous observations and modeling thus missed the hot water reservoir that dominates the new observations and inferred a much larger and colder water reservoir than we find using JWST.
Figure 9 shows the Spitzer IRS and MIRI MRS (downsampled and convolved to the IRS resolving power of 600) spectra of AS 209.We subtract the continuum baseline from the IRS spectrum using the same pipeline applied to the MIRI MRS data.The error bars in the Spitzer IRS spectrum correspond to the flux uncertainty produced by the reduction pipeline used by Pontoppidan et al. (2010b).Note that the continuum emission at both epochs is quite similar, with the Spitzer IRS continuum being at most 10% brighter near 20 µm.On the other hand, we find substantial differences between the continuum-subtracted spectra.While some spectral regions, such as between 12 -14 µm are well matched, others require the MIRI MRS spectrum to be scaled up by factors of 2 -4 to match the Spitzer IRS data.For instance, the region near the brightest CO 2 branch at 15 µm appears to be 4× as bright in Spitzer IRS, while some features between 16 -17 µm are twice as bright in Spitzer IRS.
The data suggest that the strength of the line emission in the disk of AS 209 has decreased from the epoch of the Spitzer IRS observations.Ideally, we would model the old observations to determine the extent of the variability, and whether the emission from some species has varied more than others.However, we are unable to obtain good fits to the Spitzer IRS data due to its low resolution and contrast.Rather, we provide an estimate of the variability of each species using the best fit models to the MIRI MRS data.
We estimate the (emitting area) scaling factor necessary to minimize the χ 2 between the downsampled and convolved best fit hot H 2 O model to MIRI MRS and the Spitzer IRS data, using the same wavelength regions as in Section 5.1.Then, we subtract the scaled hot H 2 O model from the Spitzer IRS data, and find whether we need an additional scaling factor to match the emission from cold H 2 O, OH, CO 2 , and HCN.Table 2 lists the total scaling factors necessary to match the emission from each species, and the wavelength ranges used to estimate each.The latter are chosen based on where the models best reproduce the MIRI MRS data, while trying to avoid regions where lines of different species are strongly blended.We find that the emitting areas of hot H 2 O and OH have decreased by a factor of 2, while those of cold H 2 O and CO 2 decreased by factors of 3 and 3.5, respectively.On the other hand, HCN does not appear to have changed.Figure 10 shows the Spitzer IRS data and the best fit model to MIRI MRS after applying the transformations described above.
It is unclear what physical mechanism could lead to the estimated emission variability.An accretion outburst during the Spitzer IRS observations could explain the brighter OH lines, as it would increase the rate of water photo-dissociation.It would also increase the disk temperature, explaining the more extended emitting areas.While there are no accurate accretion luminosity measurements at the time of the Spitzer IRS observations, based on the correlation derived by Banzatti et al. (2023b), a 3× increase in L acc would be necessary to explain the 2× higher hot H 2 O line luminosities in the Spitzer IRS observations.However, a warmer disk would also produce a brighter continuum at the Spitzer IRS epoch, yet we do not see significant differences in the dust emission.At most, the continuum in the Spitzer IRS spectrum is 10% brighter than in the MIRI spectrum, and in some regions they are fully consistent.For comparison, the continuum emission of the well-studied EX Lup disk increased by a factor of ∼4 during an outburst which increased L acc by 4.5 (Banzatti et al. 2012).
An alternative scenario could involve a 'leak' or sudden release of a large amount of icy material from an inner disk dust trap at the time of the Spitzer IRS observations (e.g.Stammler et al. 2023).The sublimation of O-rich ices could then explain why the line emission from cold H 2 O and O-carriers were preferentially affected.Otherwise, the sudden release of volatiles may be related to a giant impact, perhaps due to ongoing planet formation in the inner disk, but this scenario is highly speculative.
Ultimately, a thorough exploration of these and other ideas will be necessary to explain the observed variability, and follow-up monitoring of AS 209 with MIRI MRS may be necessary to the accomplish this goal.Furthermore, we recommend that similar comparisons with Spitzer IRS should be carried out as other sources are observed with MIRI MRS (e.g.Schwarz et al. 2023), which can help determine if this kind of line emission variability is a common phenomenon.

SUMMARY
We present new observations of the AS 209 disk using JWST MIRI MRS as part of a mid-IR follow-up to the ALMA MAPS Large Program.Our main findings are as follows: 1.The mid-IR spectrum of the AS 209 disk exhibits bright emission from hot (700 -900 K) water vapor and OH.CO 2 , HCN, and a colder water vapor component are marginally detected, whereas C 2 H 2 is not detected.We emphasize that the spectrum is severely affected by noise artifacts that prevent an accurate characterization of the molecular gas.
2. Compared to high quality MIRI data of the T-Tauri disks CI Tau and GK Tau, the hot water vapor and OH emission spectra in AS 209 do not have atypical properties.HCN and CO 2 emission features are not particularly weaker in the AS 209 spectrum.By contrast, C 2 H 2 line emission in AS 209 is weaker.
3. A comparison with Spitzer IRS observations reveals molecular emission variability in the inner disk of AS 209.The emission from hot water vapor and OH decreased by a factor of 2, cold water vapor by a factor of 3, while CO 2 emission decreased by a factor of 3.5.HCN does not appear to be affected.We consider an accretion outburst to be unlikely, given the large change in line fluxes but negligible variation in continuum emission.An explanation for the observed variability re-mains undetermined, but could be related to a triggered sublimation event and sudden release of water and organics during the epoch of the Spitzer observations.Even after the removal of strong high-frequency fringes using asteroid spectra, we find it unfeasible to visually differentiate true line emission from artifacts, and thus to subtract the continuum emission in the AS 209 spectrum using splines as is customary.We thus develop an iterative Gaussian Processes (GP) continuum-subtraction technique, which does not require the prior manual identification of line-free regions.The procedure is described below.
First, we mask out negative intensity spikes in the spectrum.To do this we use the peak finding algorithm included in scipy (scipy.signal.findpeaks).Based on visual inspection of the results, we find good performance by setting a prominence parameter of 0.05 and distance between peaks of 2 pixels2 based on visual inspection of the data.Second, we condition a Gaussian Process to the masked data.In particular, we use a linear combination of three squared exponential covariance functions (i.e.kernels) k 1 , k 2 , k 3 , such that each k n is defined as: where α and λ are the amplitude and correlation length scale, respectively, and (x i , x j ) represent any pair of points in the wavelength grid.With k 1 , we aim to capture the overall shape of the flux baseline, and we set α 1 = 100 mJy and λ 1 = 1.0 µm.Then, k 2 is set to model dust features and artifacts with α 2 = 10 mJy, but we constrain the length scale λ 2 = 0.5 µm to prevent it from erroneously removing broad emission features produced by the overlap of lines.k 3 is finally used to model the molecular line emission by setting α 3 = 1 mJy and λ 3 = 0.001 µm.These values are chosen based on prior predictive checks and visual inspection of the decomposition result.We experiment using different kernels like Matern covariance functions; the resulting decomposition is not overly sensitive to the exact choice of kernels as long as we use a mixture of functions with long and short correlation length scales, and ultimately we choose the exponential squared for its simplicity.The model is conditioned on the data using the GP Python package tinyGP (Foreman-Mackey et al. 2022) assuming, initially, a mean of zero.By setting the mean of the process to zero, the conditioned continuum function described by k 1 + k 2 does not trace a baseline under the line emission, but rather it follows the mean of the spectrum.Consequently, the line emission modeled by k 3 is incorrectly centered around zero.To solve this, we set all points of the conditioned k 3 contributions with negative values to zero, and repeat the GP conditioning using the contributions from k 3 as the new mean function of the process.At each iteration, the location of the continuum baseline is thus shifted downward by a small amount, but since λ 1 and λ 2 are constrained, the baseline is prevented from over-fitting the data and shifting the entire continuum-subtracted spectrum above zero.The process is repeated until the change in the standard deviation of the negative part of the conditioned k 3 model falls below 0.5 mJy (5 iterations).Finally, the resulting continuum baseline is linearly interpolated back onto the full wavelength grid and subtracted from the MIRI MRS spectrum, including those pixels flagged in the first step.We emphasize that the conditioned k 3 itself is not used as the continuum-subtracted spectrum, but rather the conditioned k 1 + k 2 kernels subtracted from the data.

Figure 1 .
Figure 1.Top: The complete MIRI MRS spectrum of the disk of AS 209.Bottom: The continuum-subtracted spectrum, following the procedure described in Appendix A. The shaded regions in each panel show the wavelength ranges that were excluded from our analysis.

17. 5
µm range, and only H 2 O and OH in the 20 -25.5 µm range.Note that the contributions of H 2 O and OH in the 12.3 -17.5 and 20 -25.5 µm ranges are fit together.For the ro-vibrational region, on the other hand, we only consider the contributions from H 2 O based on the rotational region results (see Section 4.1

Figure 2 .
Figure 2. Top: A closer view of the continuum-subtracted spectrum in the ro-vibrational (4.9 -8.5 µm) region.A verticallyshifted template water vapor model manually scaled to match the data is shown in violet.Middle and bottom: Same as top panel, for the rotational region of the spectrum.Note that the template model had to be scaled down by a larger amount in the ro-vibrational region compared to the rotational region.Emission lines from other species that could be visually identified are labeled in gray.

Figure 3 .
Figure 3. Bayesian Information Criterion (BIC) obtained by adding each species to the rotational region flux model.The values have been normalized for clarity.

Figure 4 .
Figure 4. Top four panels: The best fit model (green line) compared to the MIRI MRS spectrum (black line) in the 12.3 -17.5 µm range.Individual contributions from each species are shown as shaded areas.Bottom: A closer view of the brightest C2H2, HCN, and CO2 branches.Note that the best-fit C2H2, HCN, and CO2 models shown are interpreted as upper limits.
Note-From left to right: Spectral region, species name, detection status, excitation temperature, column density, emitting area, observable mass, and mass ratio with respect to hot water vapor.Best-fit parameters correspond to the median of the posterior distributions.Confidence intervals correspond to the 0.0015 and 0.9985 percentiles.Brackets indicate parameters that were kept fixed during the fit.Observable masses are the product of the best-fit column density and emitting area, times the molecular weight of each species.We do not present a mass estimate for ro-vibrational H 2 O, since the emission is affected by non-LTE excitation effects.

Figure 6 .
Figure 6.The best fit water vapor model (blue) compared to the data (black) in the ro-vibrational region.

Figure 7 .
Figure 7.The scaling factors (circles) necessary to match the optically-thick high-energy H2O line fluxes in GK Tau and CI Tau with those of AS 209.The squares and error bars show expected values and their uncertainty based on the Lacc of each host.

Figure 8 .
Figure 8.Comparison between the AS 209 spectrum and those of CI Tau and GK Tau, after scaling the latter to match the optically thick water vapor lines near 16 µm.The AS 209 spectrum is shown twice for each comparison and is shifted vertically for clarity.

Figure 9 .
Figure 9. Top: The MIRI MRS and Spitzer IRS spectra of AS 209.Bottom two panels: Continuum-subtracted Spitzer IRS and MIRI MRS spectra, convolved down to R ≈ 600.The expected location of the brightest hot (800 K) and cold (400 K) water vapor lines, as well as the brightest CO2 and organic branches are indicated.

Figure 10 .
Figure 10.The Spitzer IRS spectrum of AS 209 (black line), compared to the best fit model to the MIRI MRS data after scaling the emitting area of each species by a constant factor.The individual contributions from each species are shown as shaded regions.

Figure 11 .
Figure 11.The MIRI MRS spectrum of AS 209 in black, and the contributions from the conditioned k1 +k2 squared exponential kernels (i.e. the continuum baseline) in red.

Figure 12 .
Figure 12.The contributions from the k3 kernel after the GP model has been conditioned to the data.The model in orange shows the first iteration, and the model in teal shows the last iteration.Note that the conditioned k3 model is shifted towards positive values.

Table 1 .
Best-Fit Model Details

Table 2 .
Variability between Spitzer IRS and MIRI MRSNote-Scaling factors necessary to match the emission from each species in the best fit model to MIRI MRS, with the Spitzer IRS data.The errors correspond to 3σ confidence intervals from χ 2 curves.