CLIcK: a Continuum and Line fItting Kit for circumstellar disks

TL;DR

CLIcK is a fast, flexible tool for simultaneously fitting continuum and line emission in circumstellar disks, enabling efficient analysis of large spectral datasets, especially for upcoming JWST observations.

Contribution

It introduces CLIcK, a novel tool combining continuum and line fitting using a simplified disk model, streamlining analysis of circumstellar disk spectra.

Findings

Successfully fitted water spectra of AA Tau disk

Retrieved water vapor properties consistent with literature

Accurately recovered input parameters in simulations

Abstract

Infrared spectroscopy with medium to high spectral resolution is essential to characterize the gas content of circumstellar disks. Unfortunately, conducting continuum and line radiative transfer of thermochemical disk models is too time-consuming to carry out large parameter studies. Simpler approaches using a slab model to fit continuum-subtracted spectra require the identification of either the global or local continuum. Continuum subtraction, particularly when covering a broad wavelength range, is challenging but critical in rich molecular spectra as hot (several hundreds K) molecular emission lines can also produce a pseudo continuum. In this work, we present CLIcK, a flexible tool to simultaneously fit the continuum and line emission. The DDN01 continuum model (Dullemond et al. 2001) and a plane-parallel slab of gas in local thermodynamic equilibrium are adopted to simulate the continuum and line emission respectively, both of them are fast enough for homogeneous studies of large disk samples. We applied CLIcK to fit the observed water spectrum of the AA Tau disk and obtained water vapor properties that are consistent with literature results. We also demonstrate that CLIcK properly retrieves the input parameters used to simulate the water spectrum of a circumstellar disk. CLIcK will be a versatile tool for the interpretation of future James Webb Space Telescope spectra.

Tables

Table 1. Table 1: Modeling results obtained when fitting the 12 −  16 ​ μ ​ m 12 16 𝜇 m 12\,{-}\,16\,\mu{\rm m} H 2 ​ O subscript H 2 O \rm{H_{2}O} lines in the AA Tau IRS spectrum.

Reference	$N_{H_{2}O}\left({\mathrm{cm}}^{-2}\right)$	T (K)	$R_e\left(\mathrm{AU}\right)$	${\chi }^2$
Carr & Najita (2008)	$6.5\pm 2.4\times {10}^{17}$	$575\pm 50$	$2.1\pm 0.1$	—
Carr & Najita (2011)	$7.8\pm 2\times {10}^{17}$	$575\pm 50$	$0.85\pm 0.12$	—
This work	${3.7}_{-2.5}^{+1.0}\times {10}^{17}$	${577}_{-45}^{+35}$	${1.1}_{-0.1}^{+0.6}$	64

Table 2. Table 2: Modeling results obtained when fitting the 19 −  32 ​ μ ​ m 19 32 𝜇 m 19\,{-}\,32\,\mu{\rm m} H 2 ​ O subscript H 2 O \rm{H_{2}O} lines in the AA Tau IRS spectrum.

Reference	$N_{H_{2}O}\left({\mathrm{cm}}^{-2}\right)$	T (K)	$R_e\left(\mathrm{AU}\right)$	${\chi }^2$
This work	${4.1}_{-0.85}^{+3.21}\times {10}^{17}$	${397}_{-25}^{+15}$	${1.5}_{-0.1}^{+0.1}$	500

Table 3. Table 3: Modeling results obtained when fitting almost the entire IRS water spectrum ( 12 −  32 ​ μ ​ m 12 32 𝜇 m 12\,{-}\,32\,\mu{\rm m} ) of AA Tau.

Reference	$N_{H_{2}O}\left({\mathrm{cm}}^{-2}\right)$	$T$ (K)	$R_e\left(\mathrm{AU}\right)$	$P$	${\chi }^2$
Salyk et al. (2011)	$6.3\times {10}^{18}$	400	0.74	0	—
This work	${1.0}_{-0.2}^{+0.2}\times {10}^{18}$	${442}_{-15}^{+15}$	${1.23}_{-0.1}^{+0.1}$	0	745
(isothermal slab)
This work	${1.27}_{-0.37}^{+0.15}\times {10}^{18}$	${368}_{-5}^{+10}$ (a)	${1.79}_{-0.1}^{+0.1}$	$-{0.43}_{-0.05}^{+0.05}$	700
(Power-law profile for $T$)

Table 4. Table 4: Retrieved water line parameters for the reference model.

$N_{\mathrm{H2O}}\left({\mathrm{cm}}^{-2}\right)$	$T_{1\mathrm{AU}}$ (K)	$R_e\left(\mathrm{AU}\right)$	$P$
${7.2}_{-0.8}^{+0.4}\times {10}^{17}$	${370}_{-8}^{+2}$	$5_{-0.6}^{+0.3}$	$-{0.28}_{-0.00}^{+0.05}$

Table 5. Table 5: Parameters of the reference model.

$M_{\star }\left[M_{\odot }\right]$	0.5
$T_{\star }\left[K\right]$	4000
$L_{\star }\left[L_{\odot }\right]$	0.92
$R_{\mathrm{in}}\left[\mathrm{au}\right]$	0.1
$R_c\left[\mathrm{au}\right]$	100
$R_{\mathrm{out}}\left[\mathrm{au}\right]$	200
$\gamma$	1.0
$\beta$	1.2
$h_c\left[\mathrm{au}\right]$	10
$M_{\mathrm{dust},\mathrm{tot}}\left[M_{\odot }\right]$	$1\times {\mathrm{\hspace{0.17em}10}}^{-5}$
$\Lambda$	0.2
$f$	0.8
$D\left[\mathrm{pc}\right]$	140

Equations

$$h=h_c\times {\left(\frac{R}{R_c}\right)}^{\beta }.$$

$$\Sigma \propto {\left(\frac{R}{R_c}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{R_c}\right)}^{2-\gamma }\right],$$

$${\rho }_{\mathrm{small}}\propto \frac{(1-f)\Sigma }{h}\exp \left[-\frac{1}{2}{\left(\frac{z}{h}\right)}^2\right],$$

$${\rho }_{\mathrm{large}}\propto \frac{f\Sigma }{\Lambda h}\exp \left[-\frac{1}{2}{\left(\frac{z}{\Lambda h}\right)}^2\right],$$

Full text

1 Introduction

Circumstellar disks are a natural outcome of the star formation process (e.g., Shu et al. 1987) and provide the raw material to form planets (e.g., Lissauer & Stevenson 2007). Hence, clarifying how disks evolve and disperse is crucial to our understanding of planet formation. Although gas dominates the initial disk mass and dust accounts for only $\sim$1%, most of our understanding of disk evolution and dispersal is based on the easier to detect dust component (e.g., Williams & Cieza 2011; Ercolano & Pascucci 2017).

Infrared (IR) spectroscopy with the infrared spectrograph (IRS, Houck et al. 2004, $R\,{\sim}\,600$) on board the Spitzer Space Telescope opened a new window to investigate the gas content of disks. It revealed that the surface of disks around young solar analogues is rich in molecular emission, including water (e.g., Pontoppidan et al. 2014), as well as atomic and ionic species (e.g., Lahuis et al. 2007; Pascucci et al. 2007). High resolution ($R\,{\sim}\,25,000\,{-}\,95,000$) ground-based follow-up of the very brightest disks demonstrated that these IR lines trace the terrestrial and giant planet forming region out to ${\sim}\,10$ AU (Pascucci et al. 2011; Mandell et al. 2012). Hence, they are crucial to map the evolution of volatiles 111Molecules and atoms with low sublimation temperatures ($<$ a few 100K), including water. while planets are assemblying.

Several interesting trends have emerged from the analysis of medium-resolution Spitzer spectra. Unlike young solar analogues, disks around very low-mass stars and sub-stellar objects have spectra with stronger emission from C-bearing molecules than water (Pascucci et al. 2009, 2013) while disks around more massive Herbig Ae/Be stars lack detectable molecular emission (Pontoppidan et al. 2010) 222The only exception is the CO2 detection toward HD 101412.. In addition, the water line fluxes between $2.9\,{-}\,33$ $\mu$m correlate with the size of the inner disk gap traced by CO ro-vibrational lines suggesting that the infrared water spectrum can be used to probe inside-out water depletion in the planet-forming zone (Banzatti et al. 2017). Because line strengths depend not only on abundance but also on excitation, modeling is necessary to retrieve physical parameters such as column densities.

A few groups have carried out in depth physical and/or chemical models of a handful of well-known disks to match the dust as well as the atomic and molecular emission (e.g., Gorti et al. 2011; Kamp et al. 2013; Blevins et al. 2016). The most recent work in this area uses a new efficient and fast line raytracer FLiTs to post-process the thermochemical models of ProDiMo (Woitke et al. 2018). Such studies, although extremely valuable to constrain the properties of specific disks, are too time-consuming to investigate typical disk properties in a star-forming region and explore how those properties evolve with time. As such, other groups have taken the much simpler approach of: first, subtracting the broad continuum emission; and then, finding the best-fit (with the minimum $\chi^{2}$) synthetic spectrum assuming a simple plane-parallel slab model of gas with a single temperature and column density (e.g., Carr & Najita 2008; Banzatti et al. 2012). Interestingly, even adopting this simpler approach, molecular abundances reported in the literature differ from a few up to an order of magnitude, see e.g., the water column density inferred for the AA Tau disk from Carr & Najita (2011) and Salyk et al. (2011). These differences result from slightly different continuum subtractions, model assumptions (e.g., local line width), and choice of lines to fit. The continuum subtraction is especially critical in rich molecular spectra as hot (several hundreds K) molecular emission lines can also produce a pseudo continuum that cannot be accounted for with current approaches.

In this work, we develop a different framework to simultaneously fit the continuum and molecular emission spectra, thus accounting for any pseudo continuum produced by hot molecular lines. Our Continuum and Line fItting Kit (CLIcK) is presented in Sect. 2. In Sect. 3, we use our tool to fit the rich Spitzer water spectrum of the AA Tau disk. We also check the quality of the fit and compare the derived gas properties with literature values that are obtained from different fitting approaches. In Sect. 4, we derive synthetic spectra of a known reference model, and fit the water lines to further demonstrate the applicability of our kit. The paper ends with a brief summary and outlook in the James Webb Space Telescope (JWST) era (Sect. 5).

2 Framework of our Continuum and Line fItting Kit

There are already several workflows to simulate dust continuum and line emission in the IR regime. For instance, Pontoppidan et al. (2009) developed a line raytracer RADLite that is well coupled with the continuum radiative transfer package RADMC (Dullemond & Dominik 2004a). Dust continuum can be obtained along with the process of line raytracing. Since gas heating and cooling balance is not considered and there is no chemical network, both the gas temperature and abundance distribution of molecules must be provided either via parametric forms or from sophisticated models. Generally, one to two hours are needed to simulate a full spectrum of ${\sim}\,1,000$ water lines in the IR wavelength domain (Pontoppidan et al. 2009).

Recently, Woitke et al. (2018) ran full two-dimensional ProDiMo thermochemical disk models to calculate the dust and gas temperature structures, dust continuum and molecular concentrations. They post-processed the ProDiMo outputs with the fast line raytracer FLiTs 333http://www.michielmin.nl/codes/flits/ to simulate the molecular line spectra in the IR. A full treatment of the chemical network requires global iterations of the disk structure, which is very time-consuming even running parallel on clusters (Woitke et al. 2016).

The huge computational cost in both workflows does not favor exploring a large range of parameters when fitting spectra. Therefore, to simulate molecular line emission, we adopt a simpler approach of a slab model in local thermodynamic equilibrium (LTE) as widely adopted in previous analysis of Spitzer spectra of circumstellar disks (e.g., Carr & Najita 2008; Salyk et al. 2011; Banzatti et al. 2012; Pascucci et al. 2013). Particularly, our slab model is built upon the one developed by Pascucci et al. (2013), with molecular parameters updated to the 2016 edition of the HITRAN database (Gordon et al. 2017). The slab of gas is characterized by the temperature $T$, the molecular line-of-sight column density $N$, and the projected emitting area $A$. We define $R_{e}$ as the radius of a circular region with this emitting area. Given the parameters {$T$, $N$, $R_{e}$}, a synthetic spectrum can be computed as detailed in the Appendix of Banzatti et al. (2012). The local line width is assumed to be only due to thermal broadening. In addition to an isothermal (single temperature) slab, we also introduce a range of temperatures through a power-law temperature profile.

As mentioned above, continuum radiative transfer simulations are commonly carried out to reproduce spectral energy distributions (SEDs). The computation time depends strongly on the optical depth: for a typical circumstellar disk it is about 20 minutes without solving for vertical hydrostatic equilibrium (Dullemond & Dominik 2004a). This time is too long to generate a large number of models, which are necessary to fit the observed continuum especially when covering a broad wavelength range. Chiang & Goldreich (1997) (hereafter CG97) derived hydrostatic, radiative equilibrium models for passive circumstellar disks. CG97 models are characterized by an optically thin surface layer of superheated dust grains and a cool interior region. While the surface layer contributes most to the IR flux, the millimeter emission is dominated by the disk interior. Dullemond, Dominik and Natta (Dullemond et al. 2001, hereafter DDN01) modified the CG97 model by including a puffed-up inner rim that produces a conspicuous bump at $2\,{-}\,3\,\mu{\rm m}$ in the SED. This rim emission was proposed to explain the near-IR bump observed in the SEDs of Herbig Ae stars.

CLIcK calculates the dust continuum based on the DDN01 model which completes in just a few seconds. Details about the numerical solution are presented in Dullemond et al. (2001). In the fitting procedure, we explore the following parameters for the continuum: stellar luminosity, disk inner radius, the vertical height of the inner rim, the power-law index of the dust surface density profile, the total dust mass, and the disk inclination. The dust grains are composed of amorphous silicate and carbon, with complex refractive indices given by Dorschner et al. (1995) and Jäger et al. (1998). Mie theory is used to calculate the dust opacities. Different abundances of silicate and carbon and grain sizes are also considered, because they have a significant impact on the shape of the dust emission features such as the $10\,\mu{\rm m}$ silicate feature. The final model spectrum is a sum of the continuum and line emission. We neglect absorption of continuum radiation by molecules. Taking water as an example, a vast majority of mid-IR water lines would be optically thin if the column density of the slab is around $10^{18}\,\rm{cm^{-2}}$, i.e., a typical value found for circumstellar disks (e.g., Carr & Najita 2008; Salyk et al. 2011). Therefore, the extinction induced by molecules will not significantly change the overall shape and level of the continuum. As shown in detailed 2D models (e.g., Bruderer et al. 2015; Bosman et al. 2017), most of the line fluxes come from a region where the dust optical depth $\tau_{\rm dust}\,{<}\,1$. Therefore, absorption of line emission by the dust is also ignored. Absorption of line radiation by the gas within the slab is taken into account in the same way as detailed in the Appendix of Banzatti et al. (2012). Under these assumptions, CLIcK is flexible enough to simultaneously include multiple molecules with different sets of slab parameters into the fitting process.

The main steps of our fitting procedure, as implemented in CLIcK, are summarized in Figure 1. Since the overall flux level of a spectrum is generally characterized by the continuum, the first task (Step 1) of CLIcK is to search for a good continuum model as the starting point. To do this, we create a coarse DDN01 model grid and identify the best fit based on the $\chi^{2}$ metric. A sparse mapping of the parameter space not always provides an acceptable fit. Therefore, the second module (Step 2) invokes the Markov Chain Monte Carlo (MCMC) algorithm to improve upon the best fit returned by the grid. Details about the MCMC implementation can be found in the Appendix of Liu et al. (2012) and Liu et al. (2013). Starting from the continuum model and user-chosen input parameters for the slab of gas, CLIcK optimizes the continuum and line together (Step 3) again adopting the MCMC approach 444This approach is advantageous in high dimensional modeling such as those explored here.. As the last step (Step 4), we conduct a Bayesian analysis of the slab parameters to derive their confidence intervals. CLIcK is intentionally devised to efficiently work on a personal computer that has just a few cores. The application of our hybrid fitting framework is a compromise between the fitting quality, required time, and limited computational resources (Liu et al. 2013).

3 CLIcK on the observed spectrum of AA Tau

AA Tau is a classical T Tauri star located in the Taurus star formation region at a distance of ${\sim}137\,\rm{pc}$ (Gaia Collaboration et al. 2016, 2018). The spectral type of the star is M0.6 and the stellar mass ${\sim}\,0.6\,M_{\odot}$ (Herczeg & Hillenbrand 2014). AA Tau has been the focus of intense multi-wavelength campaigns due its periodic dimming which is ascribed to an extended non-axisymmetric feature, such as a disk warp, passing in front of a star+disk system seen close to edge-on (e.g., Bouvier et al. 2013). Recent ALMA continuum images at ${\sim}\,0.2^{\prime\prime}$ resolution have revealed that the outer disk is rather inclined with respect to the observer ($59^{\circ}$) but not edge-on (Loomis et al. 2017), providing evidence for an inner vs. outer disk misalignment as responsible for the most recent long-duration dimming that started in 2011 (e.g., Rodriguez et al. 2015). Before this dimming, AA Tau was observed with the Spitzer/IRS at a spectral resolution of $R\,{\sim}\,600$. Carr & Najita (2008) were first to report its molecular spectrum from 10 out to $37\,\mu{\rm m}$ which is particularly rich in water emission lines. As such, AA Tau is an excellent target to test and showcase our continuum and line fitting kit.

3.1 Fit to the short-wavelength IRS spectrum

The short-wavelength IRS module covers from $\sim$10 out to 20 $\mu{\rm m}$. Carr & Najita (2008) and Carr & Najita (2011) modeled the water emission spectrum of AA Tau in the wavelength range $12\,{-}\,16\,\mu{\rm m}$ using a slab of gas in LTE. They first subtracted the observed spectrum with a user-defined continuum and then performed the line fitting. Our approach differs from Carr & Najita (2008) and Carr & Najita (2011) in that we fit the continuum and line emission together.

First, we ran a grid of dust continuum SED models to identify the one that roughly matches the observed flux density versus wavelength. Then, starting from the best continuum model, together with a set of initial parameters for the $\rm{H_{2}O}$ slab, we find the best fit (with the minimum $\chi^{2}$) to the continuum and water lines by calculating 8 MCMCs, each of which containing 2,000 models. The parameters determining both the continuum and water lines are allowed to vary during the optimization. Since we are more interested to interpret the line emission, data points at wavelengths dominated by emission lines from other species (for instance HCN, $\rm{C_{2}H_{2}}$, $\rm{CO_{2}}$ and OH) or by continuum alone do not contribute to the resulting $\chi^{2}$. Our best-fit model to the $12\,{-}\,16\,\mu{\rm m}$ spectrum of AA Tau is shown in black in the upper panel of Figure 2. The lower panel of Figure 2 presents the difference between the best fit model and the observed spectrum, with the gray shaded regions indicating the wavelengths that contribute to the $\chi^{2}$ calculation. Note how the model reproduces well the observation, with relative differences typically less than 5% for the full range of fitted wavelengths.

To determine the confidence intervals of the best-fit parameters, we ran a dense grid of models with the gas temperature varying between $T_{\rm best}\,{-}\,150\,\rm{K}$ and $T_{\rm best}\,{+}\,150\,\rm{K}$ (in steps of 5 K), column density between ${\rm Log}(N_{\rm best})\,{-}\,1.5$ and ${\rm Log}(N_{\rm best})\,{+}\,1.5$ (in logarithmic steps of 0.05) and emitting area between $R_{\rm e,best}\,{-}\,1.5\,\rm{AU}$ and $R_{\rm e,best}\,{+}\,1.5\,\rm{AU}$ (in steps of 0.1 AU), where the subscript “best” refers to the best-fit values. Then, a Bayesian analysis was performed to calculate the probability distribution for each parameter (Figure 3), from which we compute 68% of the area under the functions to assign 1$\sigma$ uncertainties. We place strong constraints on the water column density and temperature, as illustrated by the clear peaks in the corresponding probability distributions and by the minor discrepancies between these peaks and the best-fit values (see Figure 3). Within 2$\sigma$, our best-fit parameters are the same as those reported in Carr & Najita (2008) and Carr & Najita (2011) (see Table 5), demonstrating the applicability of our fitting kit.

3.2 Fit to the long-wavelength IRS spectrum

Using the same optimization process described above, we fit the water lines shown in the IRS spectrum from 19 to $32\,\mu{\rm m}$. The best-fit synthetic spectrum is in good agreement with the observed one, see Figure 4. We find that the best-fit temperature ($T\,{=}\,397^{+15}_{-25}\,\rm{K}$, Table 2) is lower than the one ($T\,{=}\,577^{+35}_{-45}\,\rm{K}$, see Table 5) obtained from fitting the spectrum at shorter wavelengths, while the emitting radius is slightly larger, though the same within the uncertainty. This is as expected because emission lines at shorter wavelengths predominately trace a hotter disk region that is closer to the central star.

3.3 Fit to the entire IRS spectrum

Finally, we tested our tool by fitting almost the entire water spectrum, from 12 to $32\,\mu{\rm m}$. This has been challenging for previous methods (see Sect. 1) as it is difficult to correctly identify the global continuum in rich molecular spectra obtained at relatively poor spectral resolution. Indeed, the only attempt to retrieve molecular parameters for the entire IRS spectrum was made by Salyk et al. (2011). Instead of identifying the global continuum and fitting the entire spectrum, Salyk et al. (2011) identified 65 relatively isolated water lines, subtracted local continua, and fit their peak emission assuming a slab of gas in LTE with one temperature. Here, we use our tool to fit the continuum and water lines adopting two temperature profiles: a constant temperature (isothermal slab) and a radially decreasing power-law profile.

The upper panel of Figure 5 shows a comparison between the best-fit synthetic spectrum for the isothermal slab and the observation while Table 3 reports the best-fit parameters. Our tool does a good job in reproducing the continuum and line emission as the relative difference between the synthetic and observed spectrum is smaller than 5% for most of the considered wavelengths. The isothermal slab best-fit temperature ($T\,{=}\,442^{+15}_{-15}\,\rm{K}$) and emitting area ($R_{e}\,{=}\,1.23^{+0.1}_{-0.1}\,\rm{AU}$) are intermediate between those obtained by fitting only the short-wavelength (see Table 5) and only the long-wavelength (see Table 2) portions of the spectrum. This is as expected and further demonstrates that our tool is working properly. We also note that our best-fit temperature and emitting region are larger than those reported by Salyk et al. (2011) while the water column density is lower, see Table 3. Salyk et al. (2011) commented that fitting continuum-subtracted spectra, instead of individual line peaks, results in higher temperature, which agrees with our finding, but also higher column density and smaller emitting area, which is the opposite of what we find. However, as the latter two quantities are degenerate 777These two parameters can be considered as scaling factors for the line strength as long as the emission lines are optically thin., increasing column density coupled with a decrease in emitting area, or vice versa, will produce model fits of similar quality.

The bottom panel of Figure 5 show the fitting results for a radially decreasing power-law temperature profile. Observations with a broad wavelength coverage are sensitive to a range of gas temperatures. Consequently, this model provides an even better fit to the data than the isothermal slab as shown by the lower $\chi^{2}$, see Table 3. Even if we take into account the fact that this new model has one more free parameter, the conclusion holds because the number of fitted data points is much larger than the dimensionality of the model. Figure 6 presents the gas and dust temperature of the best-fit model, both of which share a similar power-law exponent. Moreover, the gas is hotter than the dust in the surface layer, which is consistent with circumstellar disk models that properly treat heating and cooling (e.g., Najita et al. 2011; Akimkin et al. 2013; Du & Bergin 2014).

4 CLIcK on a reference model

We further test our kit by generating and then fitting a model spectrum for a circumstellar disk located in the Taurus star-forming region. The properties of the reference model are characterized by a set of parameters which we set as a prior in the simulation. We will then discuss how well our tool reproduces the spectrum and retrieves the properties of the water emitting region.

The reference model consists of a star surrounded by a disk populated with dust and gas. Details about the dust density structure, dust properties, and stellar parameters are given in Sect. A.1. We first conduct continuum radiative transfer to self-consistently calculate the temperature structure of the dust using the RADMC package (Dullemond & Dominik 2004a). Previous observations and modeling suggest that dust and gas are thermally decoupled particularly in disk surface layers (Kamp & Dullemond 2004; Akimkin et al. 2013; Blevins et al. 2016). Therefore, we scale the gas temperature relative to the dust temperature using the results from the thermo-chemical model of Najita et al. (2011), see Sect. A.2 for details. The resulting gas temperature, together with an abundance distribution of relevant molecular species, form the input for the line radiative transfer that is performed with the RADLite code (Pontoppidan et al. 2009). Following previous tests, we will only consider water lines in the simulation, see Sect. A.2 for further details on the water abundance distribution and line radiative transfer. Finally, the output spectrum from RADLite is convolved with a spectral resolution of 2,000, the same offered by the Mid-Infrared Instrument (MIRI) on board JWST.

Using the same procedure described in Sect. 3, we fit the water lines of the reference model, with the results presented in Figures 7 and 8 and the best-fit parameters summarized in Table 4. The black line in the upper panel of Figure 7 shows the best fit synthetic spectrum which indeed reproduces most of the water lines with relative differences less than 5%. Figure 9 shows the gas temperature and water vapor column density of the reference model. On that figure we also draw the surface for dust optical depth at $20\,\mu{\rm m}$ equal to one ($\tau_{20\,\mu{\rm m}}\,{=}\,1$). The shaded area encloses the contours of ${\rm Log}\,N_{\rm H_{2}O}=17$, $T=160,\,500\,\rm{K}$ and $\tau_{20\,\mu{\rm m}}\,{=}\,1$. Therefore, these regions are expected to contribute to most of the water vapor emission at mid-IR wavelengths. Our best-fit column density ($N_{\rm H_{2}O}\,{=}\,7.2\times 10^{17}\,{\rm cm^{-2}}$), emitting area ($R_{e}\,{=}\,5\,\rm{AU}$), and temperature ($T\,{=}\,370\,\rm{K}$, $P\,{=}\,-0.28$) match well with the water vapor properties in the disk regions where most of the emission comes out, see Figure 9. This agreement shows that our tool is capable of correctly retrieving the properties of the emitting gas, even though it cannot constrain the vertical distribution of the molecules. The water vapor mass in the emitting region is ${\sim}\,40\%$ of the total water vapor mass in the reference model, in agreement with the fact that IR observations only probe a portion of the disk, e.g., the disk surface layers.

5 Summary and outlook

We have developed a flexible tool (called CLIcK) to simultaneously fit the continuum and line emission from circumstellar disks. The DDN01 continuum model and a slab model of gas in LTE are used to simulate the continuum and line emission respectively, both of them are fast enough to enable large parameter studies. We applied CLIcK to fit the observed water spectrum from the AA Tau disk and found water vapor properties that are consistent with literature values obtained by fitting continuum-subtracted spectra with a slab of gas in LTE. However, CLIcK has two advantages over these past approaches: it does not require to pre-define a global or local continuum and, in the fitting procedure, it accounts for any pseudo-continuum arising from hot molecular emission lines. We have also demonstrated that CLIcK can reproduce well the entire Spitzer/IRS water emission spectrum from AA Tau, especially with the power-law temperature profile option, and correctly retrieve the parameters of the water emitting region from a reference spectrum simulated with RADLite.

The framework of CLIcK can be easily generalized to solve similar problems in which the dust continuum has to be properly treated, including fitting the emission features from crystalline silicates (e.g., Juhász et al. 2010) and polycyclic aromatic hydrocarbons (e.g., Seok & Li 2017) in circumstellar disks.

Given the different spectral analysis used in the literature, re-fitting Spitzer/IRS archival spectra of circumstellar disks with CLIcK will be very valuable to quantify the water content in a homogeneous way. JWST will bring a remarkable leap forward in infrared spectroscopy. The spectrometers on board JWST, e.g., MIRI (Rieke et al. 2015) with a spectral resolution of $R\,{\sim}\,2,000$, will not only separate many blended lines but also boost line-to-continuum ratios, resulting in a better detection of individual lines. CLIck will be a useful tool to efficiently fit the water spectrum, subtract it, and then search for additional molecules that might be present in the atmosphere of inner disks. This will help building the chemical inventory of the terrestrial planet forming region around young stars. Moreover, with two orders of magnitude increase in sensitivity, JWST will be able to significantly increase the sample of faint disks around low-mass stars and brown dwarfs that could be observed with Spitzer. A homogeneous analysis of these spectra with CLIcK will clarify any dependencies between gas properties and stellar mass down to the brown dwarf regime.

Acknowledgements.

We thank the anonymous referee for very constructive comments that improved the manuscript. YL acknowledges supports by the Natural Science Foundation of China (Grant No. 11503087) and by the Natural Science Foundation of Jiangsu Province of China (Grant No. BK20181513) and . We thank John S. Carr for providing the reduced Spitzer/IRS spectrum of AA Tau and Klaus Pontoppidan for help with the RADLite code. We acknowledge Ewine F. van Dishoeck for insightful discussions. This material is based upon work supported by the National Aeronautics and Space Administration under Agreement No. NNX15AD94G for the program “Earths in Other Solar Systems”. The results reported herein benefited from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate.

Appendix A Simulate the infrared spectrum of the reference model

To simulate our reference spectrum we used a continuum and a line radiative transfer code. In this Appendix, we provide details about both codes.

A.1 Continuum Radiative Transfer

Density structure:    we employ the RADMC package to solve the problem of continuum radiative transfer assuming a two-dimensional flared disk model (Dullemond & Dominik 2004a). The model includes two distinct grain populations that differ in terms of dust opacity and vertical extension in the disk, i.e., a small grain population and a large grain population. The disk extends from an inner radius of 0.1 au to an outer radius of 200 au. The total dust mass in the disk is defined as $M_{\rm dust,tot}$.

Small dust grains occupy a minor fraction of the dust mass, $(1-f)\,M_{\rm dust.tot}$, and are distributed vertically following a power-law profile

[TABLE]

The degree of disk flaring is quantified with $\beta$ and $h_{c}$ represents the scale height at a fixed characteristic radius of $R_{c}\,{=}\,100\,\rm{au}$. Conversely, most of the dust mass is composed of large dust grains, $f\,M_{\rm dust.tot}$, which are concentrated close to the midplane with a scale height $\Lambda\,h$. The parameter $\Lambda\,{<}\,1$, describing the degree of dust settling, is assumed to be constant over radius $R$. This is a simple approach to account for the effect of dust sedimentation in protoplanetary disks (e.g., Dullemond & Dominik 2004b; D’Alessio et al. 2006).

We model the surface density profile as a power law with an exponential taper

[TABLE]

where $\gamma$ is the gradient parameter.

The two-dimensional density structure for each grain population can be parameterized as

[TABLE]

[TABLE]

where $\rho_{\rm{small}}$ is the density for the small grain population and $\rho_{\rm{large}}$ is the density for the large grain population. Our parameterization of the density structure has been used to successfully model multi-wavelength observations of protoplanetary disks in the literature (e.g., Andrews et al. 2011; Fang et al. 2017; Fedele et al. 2018). Table 5 summarizes the input parameters and the values we chose for the reference model.

Dust properties:    the dust in the model is a mixture of 75% amorphous silicate and 25% carbon, with complex refractive indices given by Dorschner et al. (1995) and Jäger et al. (1998). Mie theory is used to calculate the dust opacities. The grain size distribution follows the standard power law ${\rm d}n(a)\propto{a^{-3.5}}{\rm d}a$ with a minimum grain size of $a_{\rm{min}}=0.01\,\mu{\rm m}$. The maximum grain size is set to $a_{\rm{max}}\,{=}1\,\mu\rm{m}$ for the small grain population and $a_{\rm{max}}\,{=}\,1\,\rm{mm}$ for the large grain population.

Stellar heating:    we assume that the disk is passively heated by stellar irradiation. The stellar parameters are adopted as the “average” properties of a T Tauri star described by Gullbring et al. (1998). The stellar temperature and luminosity are $T_{\star}\,{=}\,4000\,\rm{K}$ and $L_{\star}\,{=}\,0.92\,L_{\odot}$ respectively. The stellar spectrum is modeled as a blackbody.

A.2 Line Radiative Transfer

As we intend to retrieve the properties of water vapor, we only simulate a mid-IR spectrum for the water lines. The line radiative transfer is performed with the RADLite code (Pontoppidan et al. 2009), which is a line raytracer designed to fast render spectra and images of complex line emission from protoplanetary disks. The molecular parameters for water are taken from the 2016 edition of the HITRAN database (Gordon et al. 2017). The level populations are assumed in local thermodynamics equilibrium, set by the local gas temperature. The local line broadening is set to $0.9\,c_{s}$, where $c_{s}$ is the sound speed.

Gas temperature:    we obtained the disk dust temperature by solving the problem of continuum radiative transfer in Sect. A.1. However, we cannot assume that the gas temperature is the same as the dust temperature throughout the disk as models that properly account for heating and cooling have showed that this assumption breaks down especially in the disk upper layers where the gas is directly exposed to ultra-violet and X-ray irradiation from the central star as well as the interstellar radiation field (e.g., Kamp & Dullemond 2004; Akimkin et al. 2013; Du & Bergin 2014). Following the approach introduced by Blevins et al. (2016), specifically their Case II, we scaled the gas temperature relative to the dust temperature. This approach is based on the thermo-chemical model presented in Najita et al. (2011), see in particular their Figure 1 which provides the scaling factor as a function of disk radius and hydrogen column density. The blue contours in Figure 9 shows the resulting gas temperature of the reference model.

Water vapor abundance:    for simplicity, we assumed a constant abundance for water vapor in the disk: $X=X_{c}(100/\rm{g2d})\,\rm{H_{2}^{-1}}$ where g2d = 100 is the gas-to-dust mass ratio and $X_{c}=10^{-4}$ is the canonical abundance for water. The abundance is set to zero at low hydrogen column densities, $N_{\rm H}\,{\lesssim}\,10^{21}\,\rm{cm}^{-2}$, to account for efficient photo-destruction (Najita et al. 2011). At temperatures lower than ${\sim}\,160\,\rm{K}$, the water is expected to condense onto dust grains. Therefore, we set the abundance at these places to a tiny value. The red dashed contours in Figure 9 show the column density of water vapor measured from the disk upper layer to the midplane of the reference model.

We focused on water lines in the wavelength range of $14\,{\sim}\,20\,\mu{\rm m}$, which is covered by the Spitzer/IRS and the upcoming JWST/MIRI instruments. The red line in Figure 7 shows the spectrum of the reference model after convolving with a spectral resolution of 2,000, similar to that offered by JWST/MIRI.

LTE approximation:    in RADLite, we assumed LTE excitation to simulate the reference spectrum, which makes the subsequent parameter retrieval fair because CLIcK calculates line emission based on the LTE assumption as well. Mid-IR water lines are dominated by pure rotational transitions. The critical densities for exciting these lines (at wavelengths considered here, $14\,{\leq}\,\lambda\,{\leq 20}\,\mu{\rm m}$) are in the range of $10^{8}\,{-}\,10^{11}\,{\rm cm}^{-3}$. In disk regions that contribute most of the water line fluxes, the gas density is comparable with the critical densities (see Figure 10), and the gas is warmer than the dust because the gas temperature is set by scaling the dust temperature based on the thermo-chemical model given in Najita et al. (2011). Therefore, LTE assumption is relatively reasonable in our case. It is beyond the scope of the paper to discuss whether or not the gas is really in LTE. Nevertheless, detailed non-LTE excitation calculations have shown that there is no dramatic deviation of the water line fluxes (e.g., Meijerink et al. 2009; Woitke et al. 2018). Similar conclusions also hold for mid-IR HCN and $\rm{CO_{2}}$ lines (Bruderer et al. 2015; Bosman et al. 2017). The abundances derived by LTE and non-LTE models do not differ by a large factor. Non-LTE excitation probably has a noticeable impact on the line shape as the emitting region can be much larger than in LTE models. However, observations with a high spectral resolving power (e.g., $R\,{\sim}\,5000$ offered by E-ELT/METIS, Brandl et al. 2014) are required to distinguish the two scenarios.