MEASUREMENTS OF WATER SURFACE SNOW LINES IN CLASSICAL PROTOPLANETARY DISKS

To map the radial distribution of water vapor, we selected a small number of protoplanetary disks known to have strong emission from water vapor at mid-infrared (10–37 μm) wavelengths, as observed with Spitzer (Salyk et al. 2008; Carr & Najita 2011), and conducted a deep spectroscopic survey with Herschel-PACS of selected water lines between 54 and 180 μm. The Herschel lines extend the Spitzer spectral mapping region from ∼1 au to ≳50 au, in particular enabling a measurement of the snow line in the line-emitting layers of the disks. To maximize the chances for detecting far-infrared water lines from disk regions with low abundances, the targeted sample was selected from disks observed to have the strongest water lines in the mid-infrared, as well as with the highest mid-infrared line-to-continuum ratios (Pontoppidan et al. 2010; Salyk et al. 2011; Najita et al. 2013). Our disk sample is composed of four classical (without known inner holes or gaps) protoplanetary disks surrounding the solar-mass pre-main-sequence stars DR Tau, FZ Tau, RNO 90, and VW Cha with spectral types spanning from M0 to G5. The water vapor observations are summarized in Table 1. The properties of the young stars and their disks are summarized in Table 2.

Table 2.  Star and Dust Disk Parameters

Name	SpT	M*	R*	Teff	Ltota	d	Mdiskb	Routb	H/Rb,c	αb,d	i	NHe	References
		(M⊙)	(R⊙)	(K)	(L⊙)	(pc)	(M⊙)	(au)			(deg)	(cm−2)
DR Tau	K7	0.8	2.1	4060	1.08	140	0.04	200	0.20	0.0357	9	4.5 × 1020	(1), (2), (3), (4), (5), (6), (7), (8)
FZ Tau	M0	0.7	2.3	3918	1.12	140	0.001	120	0.15	0.0357	45	5.5 × 1021	(5), (7), (9), (10)
RNO 90	G5	1.5	2.0	5770	4.06	125	0.005	200	0.10	0.00143	37	6.8 × 1021	(4), (5), (7), (11), (12)
VW Cha	K5	0.6	3.1	4350	3.08	178	0.001	200	0.15	0.107	45	3.6 × 1021	(4), (7), (13), (14), (15)

Notes.

^aThe luminosity is the total energy rate input into the RADMC model. It includes both stellar and accretion luminosity, but is assumed to all be launched from the stellar surface. ^bFree disk parameters. ^cDisk scale height at R_out. ^dThe power-law index describing the vertical scale height H/R ∝ R^α of the disk. ^eInterstellar extinction column densities, assuming R_V = 5.5 and ${A}_{V}/{N}_{{\rm{H}}}=5.3\times {10}^{-22}\;{\mathrm{cm}}^{2}$ (Weingartner & Draine 2001).

References. (1) Brown et al. 2013; (2) Pontoppidan et al. 2011; (3) Salyk et al. 2008; (4) Salyk et al. 2011; (5) Ricci et al. 2010; (6) Valenti et al. 1993; (7) Pontoppidan et al. 2010; (8) Hartmann et al. 1998; (9) Herczeg et al. 2009; (10) Rebull et al. 2010; (11) Levreault 1988; (12) Pontoppidan & Blevins 2014; (13) Bast et al. 2011; (14) Natta et al. 2000; (15) Guenther et al. 2007.

Download table as:  ASCIITypeset image

Deep chopping/nodding spectral range scans (R = 940–5500 or ${\rm{\Delta }}v=55\mbox{--}320\;\mathrm{km}\;{{\rm{s}}}^{-1}$) were obtained with the Photoconductor Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on board the Herschel Space Observatory (Pilbratt et al. 2010) (see Figure 2). The observations include four range scans for each disk, each providing two simultaneous spectral orders for a total of eight spectral segments per disk, five of which target rotational water transitions near 60, 75, 90, 108, and 179.5 μm (see Table 1). These spectral regions were selected to include intermediate- to low-excitation water transitions with upper-level energies spanning 114–1414 K, to complement the higher-energy (1000–4000 K) water vapor lines detected in the IRS spectra (see Figure 1). Emission from two CO lines near 90.2 μm (Meeus et al. 2012) and 108.8 μm (Karska et al. 2013) and the 5/2–3/2 OH doublet centered at 119.3 μm were also observed, but not analyzed in this paper.

Zoom In Zoom Out Reset image size

We reduced the spectra using the Herschel Interactive Processing Environment (HIPE) version 13.0.0 and extracted the spectra from the central spaxel from the rebinned data cubes. The extracted spectra were corrected to account for the flux that falls outside of the central spaxel using the HIPE aperture correction procedure, pointSourceLossCorrection.

All four disks were also observed using the Herschel-SPIRE Fourier Transform Spectrometer (Griffin et al. 2010). The SPIRE spectra were obtained in order to provide far-infrared spectrophotometry for dust modeling, as well as to measure the warm CO content of the disks. The SPIRE spectra and CO lines were analyzed in detail by van der Wiel et al. (2014).

We use R ∼ 700 (${\rm{\Delta }}v=400\;\mathrm{km}\;{{\rm{s}}}^{-1}$) mid-infrared (10–37 μm) spectra obtained with the Infrared Spectrograph (IRS) on board the Spitzer Space Telescope and first presented by Pontoppidan et al. (2010) (RNO 90, DR Tau, and VW Cha) and Najita et al. (2013) (FZ Tau). The Spitzer spectra were reduced using the pipeline described in Pontoppidan et al. (2010), following the methods by Carr & Najita (2008). Briefly, the reduction uses dedicated off-source sky integrations and an optimal selection of relative spectral response functions using all available calibration observations of standard stars. This method minimizes residual fringing and prevents the defringing procedure from removing line power from the dense molecular forest. It also enables an efficient removal of cosmic-ray hits and bad pixels.

Owing to the relatively low spectral resolving power of Spitzer-IRS, the observed water emission lines are typically blended complexes of multiple transitions (Pontoppidan et al. 2010). Further, there may be systematic errors that are difficult to quantify, especially near the edges of orders. Therefore, rather than fitting our model to the full spectral range, we opted to select 25 line complexes that have minimal systematic errors and contamination from emission by other species, including, but not limited to, OH, C₂H₂, HCN, and CO₂ (Carr & Najita 2008, 2011; Banzatti et al. 2012). We experimented with fitting the full spectral range and found that the results did not change significantly. In the Herschel ranges, individual water lines are typically well separated, with clean continua.

Integrated fluxes of Spitzer line complexes, as well as isolated Herschel lines, were estimated by fitting a first-order polynomial to the local continuum, spanning a few resolution elements on each side of the lines, and integrating the flux density in the continuum-subtracted spectra using five-point Newton–Cotes integration. Figure 14 in Appendix A summarizes the spectral regions used for the analysis. We estimated the formal statistical errors on the line fluxes by propagating the channel errors through the line integrals. In addition to the statistical, photon-noise-dominated errors, σ_phot, the high signal-to-noise ratio (S/N) Spitzer spectra often have significant, even dominant, systematic sources of error. These may include residual fringes and errors in the spectral response function. To estimate the systematic errors, we selected high-S/N spectra from bright Herbig and transition disks from Pontoppidan et al. (2010) without detected molecular line emission and obtained a smoothed estimate of the continua, F_C, using a 40-pixel-wide median filter. The effective contrast, that is, the standard deviation of the ratio of the spectrum, F, with the smoothed continuum, $C=\sigma (F/{F}_{{\rm{C}}})$, was estimated as a function of spectral channel using a moving box. The systematic error in a given spectrum was then estimated as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{sys}}(F)\sim F\times \sqrt{{C}^{2}-{({\sigma }_{\mathrm{phot}}(F)/{F}_{{\rm{C}}})}^{2}}.\end{eqnarray} \tag{ 1 }$$

Adding the systematic error to the formal statistical error in quadrature then provides the total pixel-to-pixel error of the Spitzer spectra used in our model fits. We required an S/N on each extracted line of $\geqslant 5$ for a formal detection. The selected lines and their integrated fluxes are summarized in Tables 4 and 5 in Appendix A.

To constrain the dust distribution of each disk, we combined broadband photometry with the Spitzer-IRS, Herschel-PACS, and SPIRE spectra to construct spectral energy distributions (SEDs; see Figure 3). The online SED fitting tool by Robitaille et al. (2006) was used for an initial exploration of the parameter space, while the fine-tuning of the model parameters to match the SEDs was done by hand to arrive at a plausible, but probably not unique, density structure. The stellar radius, disk mass, outer radius, flaring index, and outer vertical scale height, as defined in Dullemond et al. (2001), were varied as free parameters, while the stellar effective temperature and inclination were kept fixed. The inner-disk radius is assumed to be truncated by dust sublimation at 1500 K. For comparison to our best-fit dust models, the outer-disk scale heights relevant for isothermal disks in hydrostatic equilibrium were estimated using the model stellar radii, masses, temperatures, and disk flaring indices. The best-fit model scale height for VW Cha is close to equilibrium, but exceeds that of equilibrium for FZ Tau, RNO 90, and DR Tau by factors of 1.7, 1.4, and 2.0. In particular, the mid-infrared continuum of DR Tau was difficult to fit within the formal hydrostatic limit. We note that DR Tau has the highest accretion rate of our sample, by about an order of magnitude ($\dot{M}=1.5\times {10}^{-7}{M}_{\odot }\;{\mathrm{yr}}^{-1};$ Salyk et al. 2013). High accretion rates lead to elevated midplane temperatures, in particular within ∼5 au, and the passive RADMC model in general underestimates the midplane temperatures. This, in turn, may lead to the best-fit disk structure compensating by increasing the vertical scale heights to values that are formally too large compared to the disk temperature.

Zoom In Zoom Out Reset image size

Given the observed SEDs, the dust temperature structure for each disk was calculated using the Monte Carlo code RADMC (Dullemond & Dominik 2004), which solves the dust radiative transfer problem given an axisymmetric dust opacity and a density distribution. We use the dust opacity model from Meijerink et al. (2009) (see their Figure 2 for a plot of the absorption and scattering coefficients). The dust is a mixture of 15% amorphous carbon, by mass, from Zubko et al. (1996) and 85% amorphous silicate. It is assumed to have a size distribution parameterized using the Weingartner & Draine (2001) model, but with a maximum grain size of 40 μm and a −2.5 power slope for the silicate grains. Some stellar and disk parameters were taken from the literature and fixed for individual systems, including the stellar mass, radius, effective temperature, and inclination angle (see Table 2). The inclination angles for FZ Tau and VW Cha are not available, and we adopt a value of 45° for both disks, which is consistent with inclinations constrained using velocity-resolved CO spectroscopy by Banzatti & Pontoppidan (2015). The continuum models additionally include interstellar extinction using the Weingartner & Draine (2001) dust model with total-to-selective extinction of R_V = 5.5.

The kinetic temperature of gas in protoplanetary disk surfaces is expected to be higher than that of the dust down to vertical column densities of $\sim {10}^{22}\;{\mathrm{cm}}^{-2}$, primarily owing to collisional heating by free electrons produced by a combination of ultraviolet and X-ray irradiation (Glassgold et al. 2004; Jonkheid et al. 2004; Kamp & Dullemond 2004). In deeper layers, the electron density decreases rapidly as the gas and dust density increase, leading to strong thermal coupling between gas and dust. It is not a priori known how much of the observed water line flux is formed in disk layers with decoupled gas and dust temperatures. Simple slab models of the mid-infrared water emission derive column densities of ${N}_{{{\rm{H}}}_{2}{\rm{O}}}\sim {10}^{18}\mbox{--}{10}^{19}\;{\mathrm{cm}}^{-2}$ (Carr & Najita 2011; Salyk et al. 2011), which, assuming canonical water abundances of ${10}^{-4}\;{{\rm{H}}}_{2}^{-1}$, correspond to total gas column densities of 10²²−10²³ cm⁻², suggesting that some emitting water is in the decoupled layer, while some is coupled to the dust. In recent detailed disk models, the gas temperature is calculated including full chemical networks and assuming thermal balance between heating and cooling processes (Glassgold et al. 2009; Willacy & Woods 2009; Woitke et al. 2009; Heinzeller et al. 2011; Du & Bergin 2014). However, such calculations are CPU intensive and depend on parameters that are often poorly constrained, such as chemical rates, ultraviolet and X-ray radiation fields, and dust properties.

In order to create model grids large enough to fit the water line spectra in detail, we simplify the gas temperature calculation by computing two sets of model grids: one grid (Case I) in which the gas and dust are thermally coupled (${T}_{\mathrm{dust}}={T}_{\mathrm{gas}}$), and one (Case II) in which the gas temperature is scaled relative to the dust temperature using the thermochemical model of Najita et al. (2011). For Case II, the gas temperature scaling is computed using the ratio of gas and dust temperature as a function of disk radius and hydrogen column density ${C}_{\mathrm{scale}}={T}_{\mathrm{gas}}/{T}_{\mathrm{dust}}(R,{N}_{{\rm{H}}})$ as given by the reference thermochemical model of Glassgold et al. (2009) and Najita et al. (2011). At low column densities (${N}_{{\rm{H}}}\lesssim {10}^{21}\;{\mathrm{cm}}^{-2}$), the water is efficiently destroyed and gas temperatures rapidly rise above ∼1000 K. In this regime, the water abundance is zero. The thermochemical model has star–disk parameters that are similar to those of our observed systems (∼0.9 L_⊙, 4000 K central star surrounded by a 0.005 M_⊙ disk). Further, it includes careful treatment of the inner 10 au of the disk, where most of our line emission originates. An example of the resulting gas and dust temperatures for one of our disks is shown in Figure 4.

Zoom In Zoom Out Reset image size

The continuum models determine the radiation field at each spatial location in the disk, as well as the local dust temperatures. These quantities, along with a water abundance structure, form the input to RADLite, a line ray tracer designed to render images and spectra of complex line emission from an RADMC model (Pontoppidan et al. 2009). RADLite uses an adaptive ray grid scheme to optimally sample all size scales in a protoplanetary disk, enabling the calculation of large model grids including thousands of water lines tracing both the innermost disk and the outer disk. The molecular parameters for water are taken from the HITRAN 2008 database (Rothman et al. 2009). We assume that the level populations are in local thermodynamic equilibrium (LTE), set by the local gas kinetic temperature. The local line broadening is set to 0.9c_s, where c_s is the sound speed, consistent with recent measurements of strong turbulence in protoplanetary disk surfaces (Hughes et al. 2011). A discussion of the potential effects of assuming LTE versus non-LTE conditions is presented in Section 5.1.

The objective of this study is to retrieve a parameterized water abundance structure for each disk, given all available water line fluxes. As a starting point, we use the density- and temperature-dependent water freeze-out model from Meijerink et al. (2009) and assume a constant ice+gas abundance ${X}_{\mathrm{max}}.$ While the dust temperature is the parameter that, at the high gas densities in protoplanetary disks, primarily drives the gas/solid phase balance of water and other volatiles, there is also a weak density dependence of the water vapor abundance. This basic model has a water vapor abundance that is high at temperatures T ≳ 130–150 K and low elsewhere (see Figure 6). Herschel observations have already demonstrated that even at low temperatures, the water vapor abundance is much higher than that expected for thermal desorption (Öberg et al. 2009; Hogerheijde et al. 2011). Consequently, we additionally impose a minimum vapor abundance of X_min in the outer disk to account for the effects of nonthermal desorption mechanisms.

A fundamental property of the model is that, while the midplane water abundance drops by many orders of magnitude at the classical midplane snow line near ∼1 au, the vapor abundance is high in the superheated disk surface out to ∼5–15 au, depending on the stellar luminosity. For Case I, the radius at which optically thin dust cools below ∼150 K is the disk surface snow line and is indicated by the 150 K coupled gas and dust temperature isotherm. For Case II, the abundance model is necessarily more complex. The water vapor abundance is set to zero at low column densities, N_H ≲ 10²¹ cm⁻², owing to efficient photodestruction, following the abundance structure of Najita et al. (2011). Conversely, a high abundance, even at low column densities, results in very strong water emission at temperatures in excess of 1000 K, in clear contradiction with both data and the thermochemical models. However, we find that in this case the dust temperature is below the freeze-out temperature at disk depths high enough to prevent photodestruction of water beyond 1–2 au. Since our aim is to retrieve the radii over which water is abundant, using the line data, we must allow for abundant water vapor in disk layers where the dust temperatures are as low as ${T}_{\mathrm{dust}}\sim 60$ K. This is illustrated in Figures 5 and 6, which show the 150 K dust isotherm as reference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We explore the possibility that the water vapor abundance drops at temperatures significantly higher than 150 K. That is, we assume that the presence of cold dust always ensures that little water can remain in the gas phase while allowing for scenarios in which gas warmer than 150 K may be dry. This is motivated by thermochemical models, which suggest that efficient gas-phase formation of water is triggered at relatively high temperatures (200–300 K; Woitke et al. 2009; Du & Bergin 2014), and the observational suggestion of Meijerink et al. (2009) that the surface water abundance drops near 1 au, well within the surface snow line. Consequently, we define a critical radius ${R}_{\mathrm{crit}}\lt {R}_{\mathrm{surface}-\mathrm{snowline}}$, beyond which the water abundance drops to X_min at all disk altitudes. If the water vapor abundance is driven by pure desorption from a constant reservoir of ice, R_crit would be expected to coincide with the surface snow line. Figure 6 shows the resulting parameterized abundance structure for a case where R_crit is well inside the surface snow line. For both Case I and II, the free parameters being constrained are the inner (X_max) and outer (X_min) water abundances and the critical radius (R_crit).

While the amount of hydrogen present in the disks is not measured, the SED fit provides an estimate of the amount of dust in the upper disk layers. Assuming a gas-to-dust mass ratio leads to an estimate of the water abundance relative to hydrogen. However, since the gas-to-dust ratio is unknown, we report the water abundances relative to a canonical gas-to-dust ratio of 100.

For Case II models the gas and dust temperatures are decoupled in warm molecular disk layers. In general, we find, given identical water abundance parameters and the same critical radius, that water line strengths are strongly increased across both the Spitzer-IRS and Herschel-PACS spectral range. This is illustrated in Figure 7, and suggests that the coupled case requires higher water abundances to reproduce the data.

Zoom In Zoom Out Reset image size

Grids of model Spitzer-IRS and Herschel-PACS spectra for Case I and II were created by varying the three free parameters. R_crit was sampled in 1 au steps, while X_max and X_min were sampled in steps of 1 dex. Each model spectrum was convolved to the local spectral resolving power of the appropriate observation, and the resulting line complexes were integrated using the same procedure as for the observed spectra. The goodness of fit of the model relative to the data was evaluated using a least-squares measure. Figures 8 and 9 show the χ² surface projected onto two-dimensional marginalized space for Cases I and II, respectively. To determine the best-fit parameters without being restricted to one of the grid values, we fit a polynomial to the χ² curves in one-dimensional parameter space (Figures 15 and 16 in Appendix B). The 68% confidence intervals are estimated using a Δχ² step of 3.5, as appropriate for a three-dimensional parameter space. This leads to reduced χ² values in the range 20–32 (see Table 3). It is likely that a least-squares approach falls short for accounting for the full complexity of comparing two-dimensional dust-line radiative transfer models to high-S/N infrared spectroscopy, but developing a more sophisticated method, such as a Bayesian Monte Carlo approach, is currently unfeasible, given the computational intensity of the models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Best-fitting Water Distribution Parameters

Name	Rcrit	$X{({{\rm{H}}}_{2}{\rm{O}})}_{\mathrm{max}}$	$X{({{\rm{H}}}_{2}{\rm{O}})}_{\mathrm{min}}$	$\displaystyle \frac{{\chi }_{\mathrm{min}}^{2}}{\nu }$
	(au)	(100/gtd ${{\rm{H}}}_{2}^{-1}$)	(100/gtd ${{\rm{H}}}_{2}^{-1}$)
Case I (Tgas = Tdust)
DR Tau	>7	1.0	1.0 × 10−7	28.5
FZ Tau	5.3 ± 1.0	1.0	1.0 × 10−8	29.7
RNO 90	6.5 ± 0.4	1.0 × 10−1	1.0 × 10−7	32.2
VW Cha	7.9 ± 1.2	1.0 × 10−2	1.0 × 10−8	24.8
Case II (Tgas > Tdust)
DR Tau	3.6 ± 0.2	1.0 × 10−4	1.0 × 10−9	19.8
FZ Tau	3.3 ± 0.2	1.0 × 10−3	1.0 × 10−10	29.6
RNO 90	11.1 ± 0.2	1.0 × 10−4	1.0 × 10−9	27.6
VW Cha	3.3 ± 0.2	1.0 × 10−4	1.0 × 10−9	25.8

Note. This table lists the results obtained from modeling with (bottom) and without (top) gas–dust temperature decoupling.

Download table as:  ASCIITypeset image

Figure 10 illustrates the sensitivity of three representative water lines to different abundance parameters. Indeed, to constrain the full radial abundance structure, lines spanning the full range of upper-level energies are needed. The higher-energy lines near 17 μm are sensitive to variations in the water abundance, insensitive to variations in the critical radius, and unaffected by variations in outer-disk water abundance. The lower-energy line near 179.5 μm is sensitive to the outer-disk water abundance but insensitive to the other two parameters. Intermediate-energy lines near 60 μm are sensitive to variations in the critical radius and inner water vapor abundances.

Zoom In Zoom Out Reset image size

Emission lines from protoplanetary disks generally do not trace the cold midplane of the disk, but are formed in a layer closer to the surface owing to the temperature inversion of the superheated surface layer, line and dust opacity, or a combination of all these effects (van Zadelhoff et al. 2001; Glassgold et al. 2004; Gorti & Hollenbach 2004). Since we search for a sharp decrease in water vapor abundance at a critical radius, it is important to consider which disk height corresponds to this critical radius.

The disk region traced by rotational water lines from warm gas is illustrated using three representative lines near 17, 60, and 179 μm (shown in Figure 10), spanning upper-level energies from 2500 to 114 K. The dominance of the surface layers in the line emission is illustrated in Figure 11. Here, the vertical layers responsible for emitting between 10% and 90% of the specific intensity in the line centers are shown for three representative water lines. The line-emitting areas are superimposed on contours showing the local water abundance. Also shown are the vertically integrated water column densities of 10¹⁸, 10¹⁹, and ${10}^{20}\;{\mathrm{cm}}^{-2}$ as functions of radius. These curves allow for direct comparison to the zero-dimensional slab models, which are commonly used in the literature for characterization of molecular emission. It is seen that the mid- to far-infrared water line emission indeed originates in the disk surface, and that it is not sensitive to the location of the midplane snow line. Beyond a few au, the 179.5 μm line-emitting region spans a much wider region because the disk dust is optically thin at these wavelengths. That is, line emission from both the near and far side of the disk contributes to the total line intensity.

Zoom In Zoom Out Reset image size

In Figure 12 the radial distribution of the specific line intensity is shown for the three representative water lines. The water line at 17.23 μm, with an upper-level energy of ∼2500 K, traces radii from the inner rim of the disk out to 1 au, which is typical for lines in the Spitzer range. The 60 μm line, with an upper-level energy of ∼1500 K, traces a separate region of the disk and increases in relative contribution until it is truncated by the critical radius. Since these two lines roughly bracket the range of upper-level energies of the Spitzer range, the RADLite models appear consistent with previously published single-temperature slab model fits to the Spitzer range. For instance, Salyk et al. (2011) report water emitting radii of 0.5–2.0 au, while Carr & Najita (2011) report radii of 0.8–1.5 for disks around typical solar-mass stars. The low-temperature lines, which the 179.5 μm line represents, primarily trace the disk beyond the snow line, where abundances are low. Indeed, although the outer-disk water abundance is 5–6 orders of magnitude lower than the inner abundance, there is a significant contribution from the outer disk to the 179.5 μm line. The general nondetection of this line implies very low outer-disk water abundances (Meeus et al. 2012; Riviere-Marichalar et al. 2012).

Zoom In Zoom Out Reset image size

Table 3 shows the best-fitting water abundance parameters for the Case I and Case II grids. For the model grid with coupled gas and dust temperatures, T_dust = T_gas (Case I), we find that the best-fit inner-disk water abundances are unrealistically high ($\geqslant 0.1\;\times (100/{\rm{g}}2{\rm{d}}){{\rm{H}}}_{2}^{-1};$ see Figure 7). Since these values are much higher than those of elemental oxygen, the Case I scenario would require that either the assumed gas-to-dust ratio or another implicit assumption of the model is incorrect. Increasing the gas-to-dust ratio to 10⁴–10⁵ would be sufficient to bring the implied water vapor abundance down to the canonical value of $\sim {10}^{-4}\;{{\rm{H}}}_{2}^{-1}$. However, since the dust mass is tied to data via the SED fit, increasing the gas-to-dust ratio uniformly at all radii would also increase the disk mass to implausibly high values of ∼1 M_⊙. Therefore, if Case I were true, significant increases in gas-to-dust ratio would apply only to the surface layers of the disk traced by the water line emission.

Alternatively, allowing the gas temperature to increase, following a thermochemical calculation (Case II), leads to very different water abundances. In general, the Case II best-fit models have canonical inner-disk water abundances, even for a gas-to-dust ratio of 100 ($\sim {10}^{-4}\;{{\rm{H}}}_{2}^{-1};$ see Table 3). The retrieved inner-disk water abundances for the two cases suggest either that the gas-to-dust ratio is very high or that the gas temperature is effectively decoupled. If both effects are significant, it would suggest that the water abundance is lower than canonical.

The outer fractional water vapor abundances are representative of the cooler Herschel water tracing the outer-disk surface beyond the critical radius. Similar to the inner-disk abundances, the outer-disk abundances are lower for the Case II grid than the Case I grid, by 1–2 orders of magnitude. The Case II abundances of 10⁻⁹−10⁻¹⁰ ${{\rm{H}}}_{2}^{-1}$ are roughly consistent with previous estimates of outer-disk water abundances (Hogerheijde et al. 2011).

The model grids generally identify a best-fitting critical radius, R_crit, at which the surface abundance of water vapor decreases by orders of magnitude. R_crit is located in the range 3–10 au, although there is a tendency for the Case II best-fit R_crit to be smaller than the Case I radii. This is plausible, since Case II is able to produce more warm gas for a given R_crit, driving the parameter to smaller values. However, this effect is much less apparent than the strong decrease in fractional water abundance at all radii for Case II. Since the Case II models are more representative of the current theoretical understanding of the disk thermochemistry, and since there are many independent lines of evidence suggesting gas–dust decoupling in protoplanetary disk surface, we consider the Case II grid the best estimator of the critical radius.

In Figure 13, we compare the Case II critical radius measurements to model expectations as a function of the total system luminosity. We also compare to the midplane and surface snow lines for a fiducial passive (no accretion heating) disk for both modeling cases. The passive disk snow line relation is derived by varying the luminosity of the central star of the RADMC model for each of the four disks, while keeping all other parameters constant. It can be seen that the surface snow lines are not strongly dependent on the disk model structure. This is expected, since the surface snow line is defined as the optically thin limit and therefore approximates the isotherm of a spherically symmetric model. The midplane snow line does depend on the structure of the disk, with systematically larger snow line radii found for disks with larger scale heights, since more puffy disks intercept a larger fraction of the stellar light, which in turn leads to increasing heating of the disk at each radius.

Zoom In Zoom Out Reset image size

The measured ${R}_{\mathrm{crit}}$ are significantly larger, by a factor of ∼4, than those expected for any of the passive disk midplane snow lines. The measured R_crit are also larger than the present-day solar snow line at ∼2.7 au (Bus & Binzel 2002). However, they are significantly smaller than the surface snow line as defined by the optically thin limit. This indicates the presence of warm, but dry, surface gas and dust at radii between 4 and 15 au.

Generating large grids of model water spectra is highly CPU intensive, and several major assumptions were made to keep the problem tractable. While we allow for gas–dust decoupling in the Case II grid, level populations are set to those of LTE conditions. It is known that LTE conditions likely do not hold in general (Meijerink et al. 2009). Yet, it is also known that LTE populations lead to excellent fits of synthetic to observed spectra (Carr & Najita 2008).

Part of the explanation for this may be that the effects of gas–dust temperature decoupling and non-LTE excitation counteract each other in the disk photosphere, as shown by Meijerink et al. (2009); high gas temperatures increase upper-level populations, but the high critical densities of the water transitions conversely lead to subthermal populations at higher disk altitudes where the gas–dust decoupling is the strongest. A significant difference in our revised model, based on the Najita et al. (2011) thermochemical gas temperatures, is that the water is effectively destroyed at high altitudes, where the gas–dust decoupling is the strongest, and where densities are too low to populate the upper levels. Our LTE model spectra are therefore close to the non-LTE case of Meijerink et al. (2009). Infrared pumping tends to lead to excitation temperatures that mimic that of the local continuum color temperature, which in turn is driven by the dust temperature, rather than that of the gas. Indeed, Meijerink et al. (2009) found that, while there are indications of gas–dust decoupling and non-LTE effects in the Spitzer spectra, assuming LTE and coupled gas–dust keep line fluxes within a factor of two. While it is important to consider non-LTE effects, if such calculations can be made tractable for grid-based model fits, we do not expect them to affect the general conclusions of this paper.

We found that, in the Case II models, the local dust temperature and the location of abundant water vapor are not mutually consistent; the dust temperature is too low to allow abundant water vapor at disk radii and heights where the model fit requires it to be. Further exploration of thermochemical models and increased gas-to-dust ratio may reconcile this discrepancy.

We also neglect effects of accretion heating in the midplane by using a passive disk model. This assumption may be justified by comparing to the standard 1D + 1D models of accreting protoplanetary disks by D'Alessio et al. (2005). Using these models, Dullemond et al. (2007) compare the regimes in which the disk midplane/surface temperature structures are dominated by irradiation or accretion heating for a star–disk system similar to those of our survey (see their Figure 3). Inspection of the figure indicates that for the highest accretion rate in our sample, of just over ${10}^{-7}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}$ for DR Tau, the disk surface is irradiation dominated beyond 2 au. For our other three disks, the disk surface is irradiation dominated beyond 1 au. The D'Alessio model disk surface temperatures refer to the dust temperature. In our gas/dust decoupled case, we would expect that accretion heating plays an even smaller relative role in the gas temperature since additional external heating sources are included (Najita et al. 2011) compared to the D'Alessio models.

The best-fitting water vapor distribution models suggest that either the gas-to-dust ratio in the disk surface, at least inside the surface snow line, must be very high or most of the water emission is formed in a layer in the disk where the gas temperature is higher than that of the dust. Both of these scenarios are theoretically plausible. High gas-to-dust ratios in the inner-disk surface are predicted by many models of grain growth and settling (Miyake & Nakagawa 1995; Tanaka et al. 2005; Fromang & Nelson 2006; Ciesla 2007) and have been supported by other observations. Using near-infrared CO absorption, Rettig et al. (2006) and Horne et al. (2012) found gas-to-dust enhancements over ISM values of up to a factor 10 (or 50 in the extreme case of AA Tau). Furlan et al. (2005) found indirectly from the distribution of mid-infrared colors of protoplanetary disks in Taurus that a high degree of dust settling is common, corresponding to gas-to-dust enhancement factors of 2–3 orders of magnitude. Citing similar reasons, high gas-to-dust ratios of up to 10,000 were also used in previous two-dimensional non-LTE models of Spitzer water emission from disks (Meijerink et al. 2009).

Conversely, there is theoretical consensus that gas–dust temperature decoupling is generally significant in disk surfaces at column densities below $\sim {10}^{22}\;{\mathrm{cm}}^{-2}$ (Henning & Semenov 2013). If thermochemical gas–dust decoupling is taken into account, our best-fit Case II models indicate that there is no strong evidence for gas-to-dust mass ratios significantly in excess of 100. Other analyses of protoplanetary disk observations indicate canonical or subcanonical gas-to-dust ratios. In transition disks, recent ALMA observations of CO coupled with detailed modeling have suggested low gas-to-dust ratios of ∼10 (Bruderer et al. 2014). Thermochemical modeling of far-infrared atomic and molecular disk tracers, as observed with Herschel, has also indicated low gas-to-dust ratios in some cases (Meeus et al. 2010; Keane et al. 2014). However, it is not clear to which degree these findings are related to an advanced evolutionary stage of the disks, and the classical disks included in this paper may not be directly comparable.

Finally, it should be noted that, since the sample of water-rich disks have been selected in part based on their high line-to-continuum ratios at Spitzer wavelengths, it is possible that the uniformly high gas-to-dust ratios of our sample are a selection effect. Indeed, at the nominal water abundance and gas-to-dust ratio of ${10}^{-4}\;{{\rm{H}}}^{-2}$ and 100, respectively, the water emission would have been below the detection limit of Spitzer throughout the Case I model grid (see Figure 7).

In Figure 13, it is seen that the critical radii fall at smaller radii than the formal surface snow line, as defined by the optically thin dust limit. That is, there is a dry surface region at radii where the dust temperature is above the freeze-out temperature. There are several potential explanations for this result. Thermochemical models predict that the water abundance, in surface regions where chemical and photodestruction timescales are short, is dominated by gas-phase pathways to water with significant activation barriers (Bethell & Bergin 2009; Glassgold et al. 2009; Woitke et al. 2009). Thus, where the gas temperature drops below ∼200 K, the water vapor abundance is predicted to decrease dramatically owing to this chemical effect rather than freeze-out. Our observations are broadly consistent with this scenario. Meijerink et al. (2009) also found that the surface snow line was found at small radii in at least one case, albeit based on Spitzer data only. They suggested an alternate scenario in which strong depletion of surface water could occur if water vapor is turbulently mixed to deeper, cooler layers in the disk, where it can efficiently freeze out onto dust grains. If the icy grains settle to the midplane, the water is effectively depleted from the surface via this so-called vertical cold finger effect. In this case, it is predicted that the surface critical radius should roughly equal the radius of the midplane snow line at ∼1 au. Figure 13 indicates that the critical radius is well beyond the midplane snow line and therefore does not directly support the action of an efficient vertical cold finger effect, which depletes surface oxygen over cold midplane regions.

Further observations of water in disks are needed to confirm and refine the observed distribution of water in protoplanetary disks. We have demonstrated that, for constraining the water surface snow line, the most relevant transitions tracing 100–200 K gas are found at 40–100 μm. However, the Herschel mission has ended, and there is currently no comparable observational capability in the far-infrared. SOFIA represents one option, but it is about an order of magnitude less sensitive than Herschel and cannot observe low-excitation water, making it difficult to observe the disks around solar-type stars with strong water emission. Disk around brighter, more massive stars do not show strong water vapor emission (Pontoppidan et al. 2010; Fedele et al. 2011). ALMA may be able to observe cool water vapor using the 183 GHz water line using the new Band 5 receivers currently being constructed. In the more distant future, SPICA may provide the next set of far-infrared water vapor spectroscopy in disks, exploring lower stellar masses, and disks with fainter lines to place the strong water emitters presented in this paper into a larger, perhaps more representative, context.

A complementary option for constraining water in protoplanetary disk surfaces is to image the 3 μm water ice feature in scattered light, as demonstrated by Honda et al. (2009). This is potentially a powerful diagnostic and may be able to confirm that the disk surface is dry (no vapor or ice) between R_crit and the surface snow line.