Using Ice and Dust Lines to Constrain the Surface Densities of Protoplanetary Disks

TL;DR

This paper introduces a new method using dust and ice lines to measure the surface density of protoplanetary disks, validated on TW Hya and supported by theoretical modeling.

Contribution

The novel approach links dust and ice lines to surface density, providing a new way to estimate disk properties and test disk evolution models.

Findings

Surface density profiles are consistent with HD gas measurements.

Theoretical CO ice line estimates match observational data.

Predicted diagnostic tests for disk surface density profiles.

Abstract

We present a novel method for determining the surface density of protoplanetary disks through consideration of disk 'dust lines' which indicate the observed disk radial scale at different observational wavelengths. This method relies on the assumption that the processes of particle growth and drift control the radial scale of the disk at late stages of disk evolution such that the lifetime of the disk is equal to both the drift timescale and growth timescale of the maximum particle size at a given dust line. We provide an initial proof of concept of our model through an application to the disk TW Hya and are able to estimate the disk dust-to-gas ratio, CO abundance, and accretion rate in addition to the total disk surface density. We find that our derived surface density profile and dust-to-gas ratio are consistent with the lower limits found through measurements of HD gas. The CO ice line also depends on surface density through grain adsorption rates and drift and we find that our theoretical CO ice line estimates have clear observational analogues. We further apply our model to a large parameter space of theoretical disks and find three observational diagnostics that may be used to test its validity. First we predict that the dust lines of disks other than TW Hya will be consistent with the normalized CO surface density profile shape for those disks. Second, surface density profiles that we derive from disk ice lines should match those derived from disk dust lines. Finally, we predict that disk dust and ice lines will scale oppositely, as a function of surface density, across a large sample of disks.

Figures

a) where they find a solid depletion of up to two orders of magnitude in their fiducial model for the outer regions of late stage disks. We therefore find that our estimate of dust depletion is consistent with their model results for evolved disks.

Equations

$$T(r)=T_0\times {\left(\frac{r}{r_0}\right)}^{-3/7}$$

$$T_0=L_{\star }^{2/7}{\left(\frac{1}{4{\sigma }_{SB}\pi }\right)}^{2/7}{\left(\frac{2}{7}\right)}^{1/4}{\left(\frac{k}{\mu G{M}_{\star }}\right)}^{1/7}.$$

$$\Sigma (r)={\Sigma }_c{\left(\frac{r}{r_c}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{r_c}\right)}^{2-\gamma }\right],$$

$$\dot{M}=\frac{M_{\text{disk}}}{t_{\text{disk}}}$$

$$\begin{array}[]{@{} r @{} c @{} l @{} }&t_{s}&{}=\displaystyle\begin{cases}\rho_{s}s/\rho c_{s}&s<9\lambda/4,\text{ Epstein drag},\\ 4\rho_{s}s^{2}/9\rho c_{s}\lambda&s>9\lambda/4,\text{Re}\lesssim 1\text{ Stokes drag}\end{cases}\end{array}$$

$$\dot{r}\approx -2\eta \Omega r\left(\frac{{\tau }_s}{1+{\tau }_s^2}\right)$$

$$t_{\text{drift}}=\left|\frac{r}{\dot{r}}\right|={\Omega }^{-1}{\left(\frac{r}{H}\right)}^2{\tau }_s^{-1}\approx \frac{\Sigma r}{v_0{\rho }_{s}s}$$

$$\Sigma (r)=\frac{t_{\text{disk}}{v}_0{\rho }_{s}s}{r}$$

$$\dot{m}={\rho }_d\sigma \Delta v$$

$$\Delta v=\alpha c_s^2\frac{t_L}{t_{\eta }}({\tau }_{s,1}-{\tau }_{s,2}{)}^2$$

$$F_{\text{adsorb}}\sim n_i{c}_s$$

$$F_{\text{desorb}}\sim N_{s,i}{\nu }_{\text{vib}}{\text{e}}^{-E_i/k{T}_{\text{grain}}}{f}_{s,i}$$

$$t_{\text{des}}=\frac{{\rho }_s}{3{\mu }_i{m}_{\text{H}}}\frac{s}{N_{s,i}{\nu }_{\text{vib}}{\text{e}}^{-E_i/k{T}_{\text{grain}}}},$$

Full text

1 Introduction

Extrasolar planetary systems display a large diversity in both orbital architecture and the physical characteristics of the planets. This diversity could be the result of late stage planetary collisions (e.g. Inamdar & Schlichting, 2016), the properties and evolution of the initial gas disk (e.g. Ginzburg et al., 2016), different initial planetary formation locations in the disk (e.g. Inamdar & Schlichting, 2015) or a combination of these factors. The immediate initial conditions of planet formation are encapsulated in the protoplanetary disks that surround young stars. However, many disk characteristics remain largely unconstrained. Recent telescopic advances, particularly the Atacama Large Millimeter Array (ALMA), have enabled exploration of disks with unprecedented spatial resolution. These advances have already given us many insights into the spatial structure of disks (e.g. Nomura et al., 2016; Pérez et al., 2016; Tobin et al., 2016; Andrews et al., 2016; Williams & McPartland, 2016). Here, we take advantage of this spatial resolution to propose new observational constraints on disk surface density, a property that is fundamental for understanding both disk evolution and planet formation.

Protoplanetary disk surface density cannot be measured directly because the majority of the disk mass resides in H2 which is a symmetric particle that does not readily emit. Instead, densities have been inferred using observations of disk dust, CO, and HD. The reliability of each of these tracers has recently been called into question, leaving open the possibility that disk surface densities are entirely unconstrained (Mundy et al., 1996; Andrews et al., 2009, 2010; Isella et al., 2009, 2010; Guilloteau et al., 2011; Bergin et al., 2013; Williams & Best, 2014; Schwarz et al., 2016).

The first and most commonly used tracer of the surface density is the disk dust mass which is typically derived from resolved continuum observations or spectral energy distribution fitting (Andrews & Williams, 2005; Calvet et al., 2002; Guilloteau et al., 2011). Using dust as a tracer of total disk surface density is fallible, however, as it requires an assumed dust-to-gas ratio. This ratio is not well constrained and can be altered from the ISM value of 10*-2* (Williams & Cieza, 2011) through processes such as grain growth and particle drift (Andrews et al., 2012) and can further have a non-uniform value throughout a disk.

CO gas has also been used as a tracer of the total gas present in disks (e.g. Rosenfeld et al., 2012). However, recent observations have called into question the typically assumed abundance of CO (10*-4* in warm molecular clouds; Lacy et al., 1994), suggesting that the existence of disk processes such as photodissociation may alter this value or that there may be a global depletion of gas phase carbon in disks such as TW Hya (van Zadelhoff et al., 2001; Dutrey et al., 2003; Chapillon et al., 2008; Schwarz et al., 2016). In addition to this unknown, many CO lines are optically thick and are therefore unreliable measures of mass. The use of CO observations thus further requires careful consideration of lower optical depth CO isotopologues to estimate the gas mass to within an order of magnitude (Williams & Best, 2014; Ansdell et al., 2016).

More recently, observations of the HD $J=1-0$ line have been used to probe the gas mass in disks (Bergin et al., 2013; McClure et al., 2016). HD is thought to be a good tracer of the total gas mass as the deuterium to hydrogen fraction is relatively well-known for objects near to the sun (Linsky, 1998). However, because HD emits at temperatures above 20 K, HD emission will only trace the warm gas and thus provides a lower limit on the total gas mass present in a disk.

Given the uncertainties that accompany these observational tracers, in this paper we choose to adopt an agnostic point of view in regards to surface density. We develop a novel method to derive this quantity through asserting that the dust line locations are determined by the microphysical process of particle drift. We use physics that has been studied extensively in previous disk models. Our contribution is to suggest a new interpretation of disk observations in light of this physics.

Previous theoretical work indicates that, for the outer regions of evolved disks, drift dominates the processes of growth and collisional fragmentation in determining the maximum particle size at a particular radial scale. Particle growth can be limited by either a lack of total time (disk age), by fragmentation, or by drift. For typically turbulent disks, the disk lifetime allows plenty of time for particles to grow. The particle size at a given radial scale is limited, however, as particles are removed from the outer disk due to particle drift at a smaller maximum size than could be removed by collisional fragmentation. The maximum particle size at a given radius is thus defined as the size for which the growth timescale ($t_{\text{grow}}$) and the drift timescale ($t_{\text{drift}}$) are equal (i.e. $t_{\text{drift}}=t_{\text{grow}}$). This is described both numerically and analytically in Birnstiel et al. (2012) and Birnstiel & Andrews (2014) as the “late phase” of disk evolution in which dust growth has progressed such that it is limited by the removal of larger grains via the process of radial drift for roughly sub-centimeter sized particles.

Thus, the maximum particle size at a given radial location is limited by particle drift which will remove all particles with a drift timescale less than the age of the system (i.e. $t_{\text{disk}}$). This indicates that the equilibrium particle size at a particular disk location is limited by the time in which that particle has been able to interact dynamically with the disk gas. The maximum disk radius where we would therefore expect to see emission from a given particle size is defined as the location for which the drift timescale (and thus the growth timescale) is equal to the age of the system (i.e. $t_{\text{drift}}=t_{\text{grow}}=t_{\text{disk}}$). These locations can be seen observationally as the point where emission drops off at an observed wavelength ($\lambda_{\text{obs}}$) where we assume that the observed particle size is roughly equal to $\lambda_{\text{obs}}$. Any particle larger than this size would result in a shorter drift and/or growth timescale and would therefore reach an equilibrium location at a shorter radial scale. We note that we expect the equilibrium particle size at a particular radial location to decrease with increasing system age such that we expect larger disk radial scales at the same $\lambda_{\text{obs}}$ for younger disks, limited only by the size of the gas disk itself.

In this paper we thus make the assumption that $t_{\text{drift}}=t_{\text{grow}}=t_{\text{disk}}$ to derive fundamental disk properties. We demonstrate that, if validated, this assertion allows us to recover the surface density distribution for observed disks. To test this central premise of our modeling we use two sets of observations that give empirical information about disk dust and ice lines.

The first set of observations is a collection of recent observations, using ALMA, the Jansky Very Large Array (JVLA) and the sub-millimeter array (SMA), of TW Hya that demonstrate that the disk radial scale is distinctly smaller at longer wavelengths (Menu et al., 2014; Cleeves et al., 2015; Andrews et al., 2012). These observations provide information about the distribution of dust grains throughout the disk. We describe these observations in terms of disk dust lines which refer to the disk radial scale as it corresponds to a particular wavelength. We use these dust lines to derive disk surface densities through equating $t_{\text{disk}}$ and $t_{\text{drift}}$ as described in Section 3.1.

ALMA also provides direct observations of disk ice lines, either through direct measurements of CO gas line emission (Nomura et al., 2016; Schwarz et al., 2016) or through indirect measurements of N2H+, which is only present in large abundance when CO freezes out (Qi et al., 2013, 2015). For TW Hya, the observation of N2H+ yields an ice line location of $\sim$ 30 AU (Qi et al., 2013) while the emission from the C18O $J=3-2$ line indicates an ice line of $\sim$ 10 AU (Nomura et al., 2016). The direct observations of CO emission also give insight into the surface density contributed by CO (Cleeves et al., 2015; Rosenfeld et al., 2012). The CO ice line location depends on the CO surface density through grain adsorption rates and particle drift. These observations therefore provide a constraint on the disk surface density and CO fraction.

After discussing representative parameters for our fiducial disk, TW Hya (Section 2), we explain how we derive disk surface densities from dust lines (Section 3). In Section 4 we provide additional tests of our model through considering the disk ice lines. In Section 5 we present a description of three observational diagnostics of our model and an application to a larger range of disk parameter space. If these diagnostics confirm our interpretation, this will provide a new way to observationally measure disk masses and surface density profiles. We conclude with a paper summary and a discussion of the presented observational diagnostics in Section 6.

2 Parameters for Fiducial Disk TW Hya

We adopt TW Hya as our fiducial protoplanetary disk because it is the nearest observed disk ($d=54\pm 6$ pc) that is nearly exactly face-on ($i\sim 7^{\circ}$; Qi et al., 2004) and hence boasts a wealth of observational data. TW Hya is a long-lived disk ($t_{\text{disk}}=3-10$ Myr; Barrado Y Navascués, 2006; Vacca & Sandell, 2011) that is likely an unusually massive representative of a class of evolved protoplanetary disks as the typical disk lifetime is thought to be a few million years (Mamajek, 2009). We note that disks are typically assumed to have the same ages as their host stars and estimates of stellar ages are subject to observational uncertainties. For our discussion of TW Hya we use an approximate age of 5 Myr, however, $t_{\text{disk}}$ could be treated as a tunable quantity as appropriate.

We assume a temperature structure for TW Hya that is dominated by passive stellar irradiation. We expect that this model holds for the outer disk and note that TW Hya may be irradiation dominated at all but the very inner radii (Dullemond et al. (2007) Figure 3, based on models from D’Alessio et al. (2006)). Our parameterization of the disk midplane temperature follows Chiang & Goldreich (1997) where the canonical temperature profile is:

[TABLE]

where the coefficient $T_{0}$ is a function of stellar luminosity and stellar mass, defined at $r_{0}=$ 1 AU, and is determined via:

[TABLE]

We adopt the following parameters for TW Hya: $L_{\star}=0.28L_{\odot}$, $M_{\star}=0.8M_{\odot}$, and $\mu=2.3\text{m}_{\text{H}}$ assuming a hydrogen/helium disk composition (Qi et al., 2013; Rhee et al., 2007). Using Equation 2, we derive $T_{0}\sim 82$ K. We note that our derived midplane temperature profile is in good agreement with the observationally constrained midplane temperature derived in Cleeves et al. (2015) as well as the upper limit on the midplane temperature from Schwarz et al. (2016). We vary the normalization of this temperature profile in Section 5.2 and discuss the effect that this has on our ice line derivations. We also note that a factor of 2 change in temperature normalization ($T_{0}$) leads to a change of a factor of $\sim 0.7$ in our surface density profile derived in Section 3.1.

Spatially resolved CO observations of TW Hya have been well fit by the following surface density profile:

[TABLE]

which follows from the self-similar solution to the viscous evolution equations as shown in Lynden-Bell & Pringle (1974) and Hartmann et al. (1998). This profile is a shallow power law at small radii and follows an exponential fall off at radii larger than the critical radius, $r_{c}$. Using an assumed CO abundance of $\sim 10^{-4}$ n${}_{\text{H}}$ (the standard CO fraction in warm molecular clouds; Lacy et al., 1994), Rosenfeld et al. (2012) derive best fit parameters for TW Hya of $r_{c}=30$ AU, $\gamma=1$, and $\Sigma_{c}\sim 0.5$.

As we move forward with our discussion of TW Hya we accept the best fit parameters for all values mentioned above, except for $\Sigma_{c}$, which relies on an assumed CO abundance. Instead, we treat $\Sigma_{c}$ as a free parameter. This is motivated by discrepancies between assumed and derived CO abundances in disks. We further note that there is also a potential discrepancy between the measured and derived mass accretion rates for TW Hya.

TW Hya has an average measured accretion rate of $\sim 1.5\times 10^{-9}$ M⊙ yr*-1* (Brickhouse et al., 2012). As a consistency check, we can use the surface density profile from Equation 3 to derive an approximation of the mass accretion rate using the following expression (Jones et al., 2012):

[TABLE]

where $\dot{M}$ is the mass accretion rate, $t_{\text{disk}}$ is the age of the disk, and $M_{\text{disk}}$ is the disk mass which we take to be the mass of TW Hya interior to the critical radius of 30 AU. This estimate for $\dot{M}$ is a rough approximation for protoplanetary disks under the assumption that the primary mode of disk evolution is accretion. For instance, a rate that is higher than this derived $\dot{M}$ would quickly evolve the disk past the current state and a lower $\dot{M}$ would indicate that a process other than accretion drove the disk to evolve into its current state.

Thus, if the disk age is a proxy for evolution timescale, Equation 4 gives an accretion rate of $\sim 10^{-11}$ M⊙ yr*-1* for TW Hya using the best fit parameters from Rosenfeld et al. (2012) – a value that is 2 orders of magnitude smaller than the observational value. This value is inconsistent with observations. However, it is important to note that the accretion rates onto pre-main sequences stars are likely variable or episodic in nature (e.g. Armitage et al., 2001; Salyk et al., 2013; Hein Bertelsen et al., 2016). For episodic accretion to explain this discrepancy, TW Hya would have to be currently undergoing an episode of high accretion – a result that is unlikely given the smooth, axisymmetric nature of the disk and its observed central cavity.

We further note that photoevaporation can also remove mass in the outer disk. Since this process reduces a disk’s accretion rate onto its star for a given disk mass, if important, it would make Equation 4 an upper limit for $\dot{M}$, making its agreement with observed accretion rates worse. The rate of photoevaporative mass loss for typical fluxes is $<10^{-10}$ M☉ yr*-1* (Alexander et al., 2006), which is less than TW Hya’s observed accretion rate, so it is likely subdominant. Thus, while the estimate of the mass accretion rate from Equation 4 is not necessarily conclusive, it nevertheless provides a reason to believe that mass in the disk may be higher than indicated by CO observations.

TW Hya also has an observational lower limit of total gas mass of 0.05 M☉ from HD measurements of the warm gas in the disk (Bergin et al., 2013). The mass estimate is inconsistent with the mass estimate from the CO observations and gives an accretion rate of $\sim 6\times 10^{-9}$ M⊙ yr*-1* which is more consistent with the observed rate.

We consider the discrepancy between the measurements of the CO emission, HD gas emission, and observed accretion rate to be additional motivation for treating $\Sigma_{c}$ as a free parameter.

3 Dust Lines

The observed extent of TW Hya is wavelength dependent, ranging from a radius of $r\approx 25$ AU at a wavelength of $\lambda=9$ mm to $r\approx 130$ AU at a wavelength of 1.6 $\mu$m as shown in Figure 1. In particular, observations at 0.87 mm show a disk size of approximately 60 AU (Andrews et al., 2012, 2016), at 1.3 mm the disk size is around 50 AU (Cleeves et al., 2015), and at 9 mm the disk size is approximately 25 AU (Menu et al., 2014). We take the observed disk size of approximately 130 AU at 1.6 $\mu$m to indicate the total radial extent of the disk as this distance matches the observed radial extent of the CO emission (Qi et al., 2013).

Recent observational work has found that the continuum emission at each wavelength exhibits a markedly sharp decrease over a very narrow radial range such that $\Delta r/r\lesssim 0.1$ (Andrews et al., 2012; de Gregorio-Monsalvo et al., 2013). Thus, while these are only approximate radial sizes for the TW Hya disk, they are adequate for our preliminary physical interpretation. In our further discussion of TW Hya we allow the radius error bars to vary by $\pm 10$ AU and find little change in our theoretical modeling. These error bars are slightly larger than those derived from treating the systematic uncertainty in distance to the TW Hya system alone ($d=54\pm 6$ pc or $d=51\pm 6$ pc van Leeuwen, 2007; Mamajek, 2005, which gives uncertainties that roughly range from $\pm$ 3 AU at $\lambda_{obs}=9$ mm to $\pm$ 7 AU at $\lambda_{obs}=.87$mm). We inflate the error bars to account for error in measuring the disk radial scale without modeling the disk visibilities. We will improve on these error estimates in future work.

Models of disk emission from particles with a range of sizes indicate that $\lambda_{\text{obs}}$ roughly corresponds to the primary particle size, $s$, contributing to the observed emission (Birnstiel & Andrews, 2014). We therefore set the particle radius $s=\lambda_{\text{obs}}$ for all that follows and leave a more detailed evaluation of the particle size distributions for future work. We call the cutoff distance for emission at $\lambda_{\text{obs}}=s$ the “dust line” for particles of size $s$.

In Figure 1 we note that the observed locations of dust lines as a function of particle size (equivalent to the disk size as a function of $\lambda_{\text{obs}}$) can be well fit by a power law, with the exclusion of the disk size at the shortest $\lambda_{\text{obs}}=$ 1.6 $\mu$m. Because the radial extent of the 1.6 $\mu$m grains matches that of the CO gas, they are limited by the total disk size and not by the equilibrium that governs disk sizes at longer $\lambda_{\text{obs}}$. We conclude that this point acts as a theoretical lower limit for the radius at which you would expect to find these grains (Rosenfeld et al., 2012; Debes et al., 2013). Given their present location, they have not yet had to time to grow or drift due to the long dynamical timescales for micron sized particles in the outer disk.

The presence of dust lines in TW Hya indicate that a physical process is removing particles larger than a set maximum size from the outer disk. Particle growth can be truncated by particle fragmentation, particle drift, or a lack of total dynamical time. If we consider the case of TW Hya, the age of the system ($t_{\text{disk}}$) is sufficiently large that there is plenty of time for particle growth to occur.

Previous work has shown that for the outer regions of evolved disks, such as TW Hya, particle drift is significant for sub-centimeter sized bodies such that these particles are not able to grow to a size that can be disrupted due to collisional fragmentation (see Birnstiel et al., 2012; Birnstiel & Andrews, 2014). This is physically intuitive as small particles in the outer disk drift faster as they grow larger (see Section 3.1).

This drift limited regime reaches an equilibrium such that the radial drift of particles imposes a size limit as large particles are removed faster than they can be replenished due to particle growth. Particles in the outer disk are also removed at a size that is smaller than the limit imposed by collisional fragmentation, for disks of standard turbulence, as the relative turbulent velocities between particles of these sizes is not sufficiently large for efficient destruction (Birnstiel et al., 2012). Thus, the maximum particle size at a given radial location in the disk is given by the size for which the growth timescale is equal to the drift timescale (i.e. $t_{\text{grow}}=t_{\text{drift}}$).

Smaller particles have longer drift timescales (see Section 3.1) which means that the dust lines in a disk will evolve with disk age such that we expect the radial location of a dust line at a specific $\lambda_{\text{obs}}$ to decrease with increasing disk age. This therefore indicates that the dust line of a disk is determined by the time in which particles are able to drift, which is necessarily the age of the system. We can thus expect the fall-off of emission from a given particle size to correspond to the location where a particle’s drift timescale is equal to the lifetime of the disk ($t_{\text{drift}}=t_{\text{disk}}$). Therefore, the system can be described at any given time by an equilibrium state in which $t_{\text{drift}}=t_{\text{grow}}=t_{\text{disk}}$. This new picture of disks is summarized in Figure 2.

Using the dust line observations of TW Hya in conjunction with our theoretical premises we can now consider the dominant physical processes of growth and drift and use these calculations to derive the total disk surface density for TW Hya as well as the dust-to-gas ratio in the outer disk. We do this through the use of our central modeling premise that $t_{\text{drift}}=t_{\text{grow}}=t_{\text{disk}}$.

3.1 Radial Drift

As discussed above, we assume that the disk radial scale as a function of wavelength is set by the distance that a particle of radius, $s=\lambda_{\text{obs}}$, can drift in the age of the system. We therefore set the drift timescale $t_{\text{drift}}$ equal to the lifetime of the system $t_{\text{disk}}$.

In a protoplanetary disk particle drift occurs because the gas disk orbits at a sub-Keplerian velocity due to an outward pressure gradient (Weidenschilling, 1977). The particles in the disk continue to rotate at a Keplerian velocity ($v_{k}\equiv\Omega_{k}r$) and experience a headwind from the gas. This headwind causes the particles to loose angular momentum and drift radially inwards (Weidenschilling, 1977; Takeuchi & Lin, 2002). The amount of drift that a particle experiences depends on how well-coupled the particle is to the gas, quantified by a dimensionless stopping time: $\tau_{s}\equiv\Omega_{k}t_{s}$ where $\Omega$ is the Keplerian frequency and

[TABLE]

(summarized in Chiang & Youdin (2010)). Here $\rho$ is the gas midplane density, $\rho_{s}=2$ g cm*-3* is the density of a solid particle, $s$ is the particle size, and $\lambda=\mu/\rho\sigma_{\text{coll}}$ is the gas mean free path where $\sigma_{\text{coll}}=10^{-15}$ cm2.

The radial particle drift velocity is:

[TABLE]

where $\eta\approx c_{s}^{2}/2v_{k}^{2}$, $c_{s}$ is the sound speed of the gas, and $v_{k}$ is the Keplerian velocity (see the review by Chiang & Youdin, 2010).

We can now derive an equation for drift timescale ($t_{\text{drift}}=\left|r/\dot{r}\right|$) that directly depends on the disk surface density, $r$, and $s$. We first note that the particle sizes we consider ($s<1$ cm; see Section 3) typically interact with the disk gas via the Epstein law for gas drag and have a dimensionless stopping time less than 1. Once our calculation is complete we verify that the Epstein regime applies for our reconstructed surface density profile and that $\tau_{s}<1$.

We approximate a dimensionless stopping time in the Epstein regime as $\tau_{s}=\Omega\rho_{s}s/\rho c_{s}$. We can rewrite this quantity by noting that the disk’s scale height $H=c_{s}/\Omega$ and its surface density in gas $\Sigma_{g}=H\rho$ to obtain $\tau_{s}=\rho_{s}s/\Sigma_{g}$. As found in Rosenfeld et al. (2013), increasing temperature with height can inflate the aspect ratio of an observed disk when layers several scale-heights above the midplane are probed. As most of the disk mass is contained within the scale height closest to the midplane, however, the temperature at the midplane sets the volumetric density in the disk to an order of magnitude as described above.

We reformulate the drift timescale through defining the parameter $v_{0}$ which physically corresponds to the maximum drift velocity. For a passively irradiated disk, $v_{0}\equiv\eta v_{k}=c_{s}^{2}/2v_{k}$ varies very weakly with radius. We find that $c_{s}\propto r^{-3/14}$ where $c_{s}=\sqrt{kT/\mu}$ using Equation 1 for temperature as a function of radius. For the Keplerian velocity we find that $v_{k}\propto r^{-1/2}$ using Kepler’s third law. This gives a parameter $v_{0}$ that is weakly dependent on radius: $v_{0}\propto r^{1/14}$.

We can now write the drift timescale directly in terms of the surface density, radius, particle size, and the maximum drift velocity $v_{0}$:

[TABLE]

We can now set the drift timescale equal to the age of the disk ($t_{\text{drift}}=t_{\text{disk}}$) and solve for the gas surface density at a given dust line radius. We thus derive the following equation for disk surface density as a function of radius.

[TABLE]

From this equation we see that $\Sigma(r)\propto s/r$ with $s(r)$ plotted in Figure 1. This reformulation is thus particularly powerful as it gives a direct scaling between surface density and radius with only a very weak dependence on the temperature profile of the disk via the term $v_{0}$.

We now apply Equation 8 to TW Hya again assuming that the dominant particle size $s$ is given by the wavelength of the observation ($\lambda_{\text{obs}}$). The resulting surface density points are well described by a steep power law function of approximately $r^{-4}$. Initially, this appears to be an unrealistically steep relation for the disk surface density as typical scalings for the minimum mass solar nebula are $\Sigma\propto r^{-1}$. However, as these points fall close to the observationally derived exponential fall-off range for the measurements of CO in Rosenfeld et al. (2012), the data match this profile quite well when we allow $\Sigma_{c}$ to be a floating normalization factor (see Section 2). We find that a surface density normalization of $\Sigma_{c}\approx 10^{2.5}$, approximately 3 orders of magnitude larger than $\Sigma_{c}$ from Rosenfeld et al., adequately matches our derived surface density points. The derived surface density profile with error estimates is shown in Figure 3. We find that the particle sizes we consider in the outer disk are indeed still well within the Epstein regime as was assumed in our derivation of the disk surface density profile.

We note that this surface density estimate depends on the age of the disk which may be uncertain by as much as a factor of 2 (e.g. Pecaut et al., 2012). For TW Hydra the surface density normalization, $\Sigma_{c}$, ranges from $\sim$ 200 at 3 Myr to $\sim$ 800 at 10 Myr. We note that the upper end of this age range will lead to a derived disk surface density that is Toomre-Q unstable, however, due to the order of magnitude nature of this derivation those numbers could be revised through a more detailed calculation. This trend generally holds such that the derived disk surface density increases with an increase in estimated disk age.

We can now convert our surface density profile to a derived gas mass as well as a mass accretion rate for TW Hya. We find that our derived gas mass of approximately 0.05 $M_{\odot}$ is consistent with the observational lower limit of 0.05 M⊙ as derived from HD measurements of the warm gas (Bergin et al., 2013). We remark that this agreement implies no low-density cool gas and note that there is room within the errors for there to be a comparable amount of cold and hot mass present within the disk, however, to avoid being Toomre-Q unstable the mass cannot increase by a substantial amount.

Following the discussion in Section 2, with the inclusion of rough radius error estimates, we find a mass accretion rate of $\dot{M}\sim 4\times 10^{-9}-2\times 10^{-8}$ M⊙ yr*-1*. Compared to the measured value of $\sim 1.5\times 10^{-9}$ M⊙ yr*-1* (Brickhouse et al., 2012) our derived accretion rates give a larger value more in line with measured accretion rates for younger systems than the derived disk accretion rate using CO surface density observations as in Rosenfeld et al. (2012). We note that TW Hya has an inner disk gap at 1 AU (Andrews et al., 2016). It is possible that this gap was formed by the presence of an accreting protoplanet which could be allowing only $\sim$ 10% of the accretion flux onto the star (Papaloizou & Nelson, 2005; Najita et al., 2007; Zhu, 2015; Eisner, 2015).

We also note that, given our estimated disk age, our derived surface density profile is Toomre-Q unstable for the higher region of our estimated errors. We note, however, that assuming a shorter age of the system ($\sim 3$ Myr) gives a stable disk surface density profile for all estimated errors.

3.2 Particle Growth and the Dust-to-Gas Ratio

Given that TW Hya is an evolved disk in dynamical equilibrium, we can use an approximation of the growth timescale to determine the dust-to-gas ratio in the outer disk through allowing $t_{\text{grow}}=t_{\text{disk}}$ as discussed in Section 3. We derive our growth timescale by first considering the growth rate in terms of particle size ($s$):

[TABLE]

where $\rho_{p}$ is the density of the particles which we take to be $\sim\rho_{g}f_{d}$ where $f_{d}$ is the dust-to-gas ratio and $\sigma=\pi s^{2}$ where $s$ is the size of the dominant particle at that radius $r$. In the outer disk we assume that the relative particle velocities $\Delta v$ are due to turbulence which we approximate as:

[TABLE]

(Ormel & Cuzzi, 2007) where $t_{L}$ is the overturn time of the largest eddies which we take to be the orbital period, $t_{\eta}=Re^{-1/2}t_{L}$, $\tau_{s,1}$ is the stopping time of our dominant particle size, $\tau_{s,2}$ is the stopping time of the particle size that contributes to growth, and $\alpha$ is the parameter of ignorance for turbulence in a disk which we take to be a standard value of $10^{-3}$ (e.g. Rafikov, 2017). This $\Delta v$ is valid for particles that are tightly coupled to the gas. We note that direct observations of disk turbulence (Hughes et al., 2011a; Teague et al., 2016) do not yet probe the disk midplane.

We are now able to solve for the dust-to-gas ratio in the disk through setting the growth timescale ($t_{\text{grow}}$) equal to the age of the system ($t_{\text{disk}}$) and solving for $f_{d}$. We do this while assuming that the particles grow through perfectly efficient sticking with particles that are similar in size. We make this approximation as most of the dust surface density will be in the largest grains such that these particles dominate particle growth. Given this assumption, typical values of $\Delta v$ for our TW Hya calculation are $\sim 3-10$ cm s*-1*. For TW Hya this exercise gives an average dust-to-gas ratio of $f_{d}\approx 10^{-3}$ in the outer disk where we observe the disk dust lines.

As a consistency check we can compare this derived dust-to-gas ratio with the observed dust surface density profile described in Andrews et al. (2012). We note that the total dust mass and the dust surface density profile are model dependent quantities with values that vary across the literature (Calvet et al., 2002; Andrews & Williams, 2005; Hughes et al., 2008, 2011b; Kamp et al., 2013; Hogerheijde et al., 2016). We thus choose the Andrews et al. (2012) model as it is based off of 870 $\mu$m emission which is within the range of our dust line observations.

In their model, the dust surface density profile is well described by a shallow power law $\Sigma_{d}\propto r^{-0.75}$ until the emission falls off at roughly 60 AU. This profile is normalized such that $\Sigma_{d}=0.39$ g cm*-2* at 10 AU. When we divide this dust surface density profile by our derived total surface density profile we find an average dust-to-gas ratio of $\sim 10^{-2.5}$ which is roughly consistent with our derived dust-to-gas ratio above.

We further compare our dust surface density points at the three dust line locations using our derived dust-to-gas ratios with the best fit model from Andrews et al. (2012) in Figure 4. We find a dust surface density that is systematically lower than the best-fit model from the 870 $\mu$m emission, however, without error bars on this observationally derived model it is difficult to know if our results are discrepant. We determine an extremely rough estimate of the uncertainty of this model by propagating uncertainties on the mass opacity, flux calibration ($\sim 10$ % systematic uncertainty), and distance (a few percent uncertainty). The largest source of uncertainty in the Andrews et al. (2012) model is the mass opacity. We adopt an uncertainty on this value of roughly $\sim 99$ % as the range in assumed mass opacities for dust at $\sim$ 870 $\mu$m ranges from $\sim$ 10${}^{-2}-10^{0}$ in the relevant literature (e.g. Hildebrand, 1983; Pollack et al., 1994; Henning & Stognienko, 1996; Andrews & Williams, 2005; Andrews et al., 2012). The error on the mass opacity thus dominates the uncertainty in the dust measurement. We note that our points represent the dust surface density comprised of the large particles in some set size distribution and that including the contribution of smaller grains to the total surface density should increase our answer by a few tens of percent.

We also compare to dust depletion in a time evolving disk as shown in Birnstiel & Andrews (2014) (see their Figure 4a) where they find a solid depletion of up to two orders of magnitude in their fiducial model for the outer regions of late stage disks. We therefore find that our estimate of dust depletion is consistent with their model results for evolved disks.

In conclusion, by considering particle drift we derive a total disk surface density for TW Hya that is roughly 3 orders of magnitude larger than the surface density derived from CO observations through considering particle drift. We further derive a dust-to-gas ratio of approximately $10^{-3}$ in the outer disk ($\sim 25-60$ AU).

4 Ice Lines

We can now further our discussion of the disk surface density through a careful treatment of the molecular ice lines. The three ice lines that are most frequently calculated are H2O, CO2, and CO as these species are considered to be in relatively high abundance (Öberg et al., 2011). To date, the CO ice line is the most readily observed due to its large radial distance from the central star (Qi et al., 2013). We thus primarily focus on the CO ice line in our calculations. However, when we expand our discussion to consider a larger parameter space of disk parameters, we further extend our arguments to the other volatile species as well.

While we consider a passively heated disk (see Section 2), we note that the water ice line may be impacted by accretion heating. We include some plots of the H2O ice line for reference, but these should be viewed with caution.

Our discussion focuses on the midplane ice lines as opposed to the surface ice lines frequently discussed in the literature (Blevins et al., 2016; Schwarz et al., 2016). As a result we will compare our derived ice lines to observations that probe the midplane ice line directly (via the C18O line; Qi et al., 2015; Dartois et al., 2003) or indirectly (via the N2H+ ion; Qi et al., 2013, 2015).

There are three pieces of physics that we consider in our ice line calculations: particle adsorption, desorption, and drift.

4.1 Volatile Adsorption and Desorption

The classic ice line calculation balances adsorption and desorption flux onto a grain to determine the ice line radius (Hollenbach et al., 2009; Öberg et al., 2011). We refer to this ice line as the ‘classical ice line’. Following Hollenbach et al. 2009, these two fluxes are quantified as:

[TABLE]

[TABLE]

where $n_{i}$ is the relevant gas density species, $c_{s}$ is the sound speed, $N_{s,i}\approx 10^{15}$ sites cm*-2* is the number of adsorption sites per volatile per cm2, $\nu_{\text{vib}}=1.6\times 10^{11}\sqrt{(E_{i}/\mu_{i})}$ s*-1* is the molecules vibrational frequency in the surface potential well, $E_{i}$ is the adsorption binding energy in units of Kelvin, and $f_{s,i}$ is the fraction of the surface adsorption sites that are occupied by species $i$ (which we take to be unity). Finally, we assume that $T_{\text{grain}}=T$, meaning that the dust and gas have the same temperature in the disk mid-plane.

Balancing Equations 11 and 12 allows for us to solve for the freezing temperature of a species as a function of radius for a given disk surface density profile. We then locate the classical ice line by finding the disk radius where the molecular freezing temperature is equal to the disk temperature. This self-consistent method allows us to determine how the ice line location changes as a function of both disk surface density and temperature.

4.2 The Influence of Particle Drift

Particle drift, as described in Section 3, influences the location of the ice lines (Piso et al., 2015, 2016). This is because particles that drift faster can cross the ice line before desorbing, thus potentially moving the location of the ice line inwards. The drift ice line location can be calculated by setting the desorption timescale,

[TABLE]

(where $\mu_{i}$ is the molecular weight of the desorbing species) equal to the drift timescale (the radial location of a particle divided by Equation 6) and solving for the desorption distance, $r_{\text{des}}$ (as verified by time evolving calculations in Piso et al., 2015). This is done analytically for the small stopping time approximation in the Appendix of Piso et al. (2015). We extend this calculation to the $\tau_{s}>1$ regime through making no approximations in Equation 6 in regards to $\tau_{s}$. Instead, we use the two distinct stopping time expressions given both Stokes and Epstein regimes (Equation 5). We use the SciPy subroutine fsolve (Jones et al., 2001) to solve for radius through balancing the drift timescale with the desoprtion timescale ($t_{\text{drift}}=t_{\text{des}}$) where these timescales are given by $t_{\text{drift}}=\left|r/\dot{r}\right|$, where $\dot{r}$ is given by Equation 6, and Equation 13. We consider the radial dependence of each term without approximation, which allows us to derive a completely self-consistent estimate of $r_{\text{des}}$ for any arbitrary surface density and/or temperature profile.

We compare the difference between the drift ice line and the classical ice line for our new self-consistent solver and the analytic solver from Piso et al. (2015) in Figure 5. We use the disk temperature and surface density profiles from Piso et al. (2015) and find that particles smaller than $\sim$ 10 cm ($\tau_{s}\sim 0.5$) do not drift past the CO ice line. Particles larger than $\sim$ 10 cm do experience drift, with the maximum drift reached at a $\tau_{s}\approx 1$ that then reaches a near constant value at higher stopping time. We note that, while small particles do not drift to the same radial location as large particles, accretion of disk gas has been shown to pull smaller particles in past the ice line as shown in Piso et al. (2015).

The nearly constant drift distance at particle sizes with a stopping time greater than 1 occurs due to the interplay between two different pieces of physics. The first is that the larger the particle size the longer it takes for that particle to desorb and the further it can drift past the ice line. The second is that the maximum drift velocity occurs for particles with $\tau_{s}=1$. The decrease in drift velocity past a stopping time of one is offset by the increase in particle size until the regime changes from Epstein to Stokes. The drift timescale for a particle in the Epstein with $\tau_{s}>>1$ is $t_{\text{drift}}=\left|r/\dot{r}\right|\sim\tau_{s}/\eta\Omega\sim\rho_{s}s/\rho c_{s}\eta$. The desorption timescale is simply proportional to the particle size ($t_{\text{des}}\propto s$) with a desorption constant such that $t_{\text{des}}\sim C_{\text{des}}s$. When we equate these two timescales we find that the dependence on particle size cancels out and that the drift distance is therefore independent of this quantity. We find that this approximation is roughly accurate for particles in the Epstein regime with $\tau_{s}$ close to or greater than 1.

Due to this effect, we consider the “drift ice line” to be the location for which a $\tau_{s}=1$ particle starting at the classical ice line location drifts and then desorbs for situations in which we know the disk surface density. If we do not know the disk surface density we use a large sized particle such that it desorbs at the maximum drift radius for a range of disk parameters. We numerically solve for the drift distance without making any further approximations in all further discussions of drift. The existence of a maximum drift distance found in our self-consistent solution allows us to predict the circumstances in which drift affects the ice line locations.

We find that the radial location of the classical ice line strongly impacts the amount in which drift will play a role. In particular, we find that the process of drift becomes important when the classical ice line is located where larger particles have a stopping time close to 1 as $\tau_{s}=1$ particles drift the fastest and larger particles take longer to desorb (see Figure 6). The size of the $\tau_{s}=1$ particle depends on the overall surface density profile of the disk. Drift will therefore matter the most in determining the true ice line location at the disk radial location where the particles with $\tau_{s}=1$ reach a maximum size. If the classical ice line of a particular molecule is located near this point then drift will affect the location of the ice line radius pushing it inwards. This is shown in Figure 6 through a consideration of our derived surface density profile for TW Hya and the standard profile used for the minimum mass solar nebula ($\Sigma\propto 1700\text{ g cm}^{-2}\times r^{-3/2}$) (Öberg et al., 2011; Chiang & Youdin, 2010).

4.3 Application to TW Hya

Given our derived total surface density profile for TW Hya, the last parameter needed to derive the CO ice line location is the CO abundance. We use the measured CO surface density from Schwarz et al. 2016 in conjunction with our derived total surface density profile to uncover an approximate CO abundance of 10*-7* n${}_{\text{H}}$ from $\sim$ $10-60$ AU. This value is consistent with the average upper limit of 10*-6* n${}_{\text{H}}$ found in Schwarz et al. (2016) which was derived by comparing their CO surface density to the surface density profile derived from HD observations of the warm gas (Bergin et al., 2013).

This reduction of 3 orders of magnitude from the measured abundance in the ISM (10*-4* n${}_{\text{H}}$) demonstrates that our model requires a global depletion of volatile carbon in TW Hya. Cleeves et al. (2015) infer a similar level of depletion using observations of CO and modeling by Du et al. (2015) supports this conclusion.

We can now calculate both the classical and maximum drift CO ice lines. For our calculations we use a binding energy of $E_{i}/k\sim 850$ K (Aikawa et al., 1996) 111We note that observations discussed in Schwarz et al. (2016) indicate a CO binding energy closer to $E_{i}/k\sim 960$ K for TW Hya. Using this binding energy we determine that the classical CO ice line is located at $\sim$ 23 AU and the drift ice line is located at $\sim$ 10 AU. These derived values are also reasonably consistent with the observed ice lines of TW Hya (see text).. For TW Hya, we derive a classical CO ice line of $\sim$ 30 AU and a drift ice line of $\sim$ 15 AU.

We find that our derived classical ice line is in agreement with the ice line determined observationally via the N2H+ ion of $\sim$ 30 AU (Qi et al., 2013). We further find that our derived drift ice line is in good agreement with the CO ice line as derived through the C18O $J=3-2$ line where a decrement of CO emission is observed until $\sim$ 10 AU (Nomura et al., 2016). Given the discrepancy between the ice line as determined by N2H+ and C18O and the corresponding agreement we find with our theoretical model, we note that the C18O $J=3-2$ line may be the most sensitive probe of a disk’s midplane ice line, particularly when particle drift effectively smears the ice line out over a relatively large radial scale. This could be because the presence of N2H+ requires a lack of CO gas in the disk and may thus require that all particle sizes have frozen out (i.e. the classical ice line).

We can now put together our complete model for the dust and ice lines in TW Hya as shown in Figure 7.

5 Observational Diagnostics

We now propose three different observational tests of our model framework beyond the application to TW Hya. The first test is whether or not other disks have dust lines that can be well fit by a normalized form of the physically motivated surface density profile in Equation 3 for disks where $r_{c}$ is constrained by CO observations. This is a particularly powerful test of concept if we find dust lines in other disks that are located in the power-law region of the surface density profile interior to the exponential fall-off. The second test is whether we can derive a surface density profile from the disk ice lines that is consistent with the surface density profile derived from disk dust lines. The third test is whether or not the dust and ice lines scale oppositely, as a function of surface density, across a large sample of disks.

Our second and third tests rely on an understanding of both classical (i.e. classical regime) and drift (i.e. drift dominated regime) ice lines as described in Section 4. The particular assumptions and/or prior knowledge of a disk will inform which of the two regimes is relevant as discussed in Section 5.2.

5.1 Test 1: Surface Density from Disk Dust Lines

The first test of our model is whether or not other disks have dust lines that can be well fit by a normalized form of the physically motivated surface density profile in Equation 3. In the era of ALMA, multiple disks will have high resolution data at several wavelengths – the ideal data set for this fundamental test. Following the technique described in Section 3 we can convert these dust lines into a surface density measurement and see if these measurements scale with the observationally derived surface density profiles.

In the case of TW Hya, the disk dust lines were near the exponential fall-off point, which unexpectedly (due to the steep slope of the derived power law) matched previous surface density models with the application of a normalization factor. For other disks there is no reason to think that the disk dust lines will necessarily fall in the same exponential fall-off region of the disk surface density profile. If a disk has a critical radius further from the star then it is increasingly likely that we will be able to detect dust lines interior to this point.

This test is thus a powerful test of concept that will allow us to see whether other disk sizes can be well fit through considering drift as the primary driver of disk dust lines and using these derived surface density values to normalize previous observationally fit disk surface density profiles.

• •

Test 1 Observational Requirements: This test can be best carried out on disks with CO observations, such that there is a derived critical radius, and observations of dust emission at several wavelengths such that the disk dust lines are known.

5.2 Ice Line Regimes

It is now helpful for us to distinguish between the drift and classical ice lines from observations of CO or other ice line tracers. Here we provide several rules of thumb to aid in the use of ice lines as useful diagnostics. However, we note that these aids should be used initially and then verified for self-consistent results.

Our proposed tests using these ice lines differ as the classical ice line depends on the CO abundance while the drift ice line does not. This makes the drift ice line particularly useful as it relies on fewer assumptions while the classical ice line retains its usefulness when considered across a sample of disks.

As mentioned in Section 5, there are two regimes that dictate when the drift or classical ice lines are relevant: the drift-dominated regime and the classical regime. In the drift-dominated regime the drift ice line is interior to the classical ice line and thus there should be only be CO freeze out exterior to this point with complete freeze-out occurring exterior to the classical ice line, making the drift ice line the most interior detectable ice line in the disk. In the classical regime particle drift does not happen quickly enough for particles to cross the ice line without desorbing and the classical ice line is observed.

To interpret the surface density from the disk ice line (see Section 5.3) we need a priori knowledge of the relevant ice line regime in the disk. Without knowledge of the disk surface density profile we can roughly determine the correct regime through an analysis of the ice line dependence on elemental abundance and disk temperature as well as a comparison of the ice line location with respect to the disk’s critical radius.

One simple metric in determining the ice line regime is the knowledge that the ice line location is preferentially dominated by drift for disks that have large particles with a stopping time of unity near the ice line location. As seen in Section 4.2, particles can drift further without desorbing where the $\tau_{s}\approx 1$ particle size is large. This size increases at a given radius with increased disk critical radii and surface density as demonstrated in Figure 8 where we calculate the $\tau_{s}=1$ particle size at 30 AU for a disk with the same temperature profile as TW Hya for a range of surface density normalization factors and disk critical radii. We find that an ice line observed interior to the critical radius is likely dominated by drift for massive disks with large critical radii.

Another parameter that determines the ice line regime is the amount of radiation that the disk receives from its host star (i.e. the stellar luminosity). The importance of this parameter is shown in Figure 9 for generally assumed molecular abundances ($\text{n}_{\text{CO}}=1.5\times 10^{-4}\;\text{n}_{\text{H}}$, $\text{n}_{\text{CO}_{2}}=0.3\times 10^{-4}\;\text{n}_{\text{H}}$, $\text{n}_{\text{H}_{2}\text{O}}=0.9\times 10^{-4}\;\text{n}_{\text{H}}$; Pontoppidan, 2006).

We find that we are in the drift dominated regime for disks with moderate densities and high temperatures across all three molecular species. At very high disk densities we find that the drift ice line is irrelevant and only the classical ice line can be observed, as a large increase in disk density moves the classical ice line inwards such that the disk temperature at the ice line is higher and the desorption timescale is shorter (i.e. the increased temperature near the star overwhelms the effect of the increase in disk density, see Equation 12). As the disk density is not known a priori we posit that the drift ice line is likely relevant for disks with high stellar irradiation.

We further find that, given a molecular depletion as in the case of our derived disk parameters for TW Hya and may be true for the majority of disks, drift sets the true ice line location across a wide range of parameter space as shown in Figure 10. Drift is more important when the CO abundance is smaller because the decrease in CO abundance moves the classical ice line further from the star while the drift ice line remains unchanged with the change in abundance. Therefore, we find that disks depleted in CO should have ice lines that are determined by drift. This is further demonstrated in Figure 11 where, across a wide range of molecular CO abundances, the drift ice line is constant and interior to the classical ice line given our derived disk parameters for TW Hya.

It is therefore a relatively safe assumption to consider the CO ice line location to be determined by drift for disks that are hot, depleted in CO, and/or have an observed ice line interior to a large critical radius and appear to be relatively massive. Disks for which the opposite is true are in the classical regime and their ice lines should be theoretically treated as classical ice lines.

5.3 Test 2: Surface Density from Disk Ice Lines

The second test of our model is whether or not the surface density profile derived from a consideration of the disk ice line matches the surface density profile derived from the disk dust lines. We now discuss how empirical knowledge of a disk’s molecular ice line can be used to independently determine the total surface density profile of a disk. This can be done as both the drift and classical ice lines depend on the disk surface density.

5.3.1 Drift Ice Line Surface Density Determination

If the disk in question is in the drift dominated ice line regime (see Section 5.2), solving for the disk surface density profile simply involves setting Equation 6 for $\tau_{s}=1$ (where particle drift reaches a maximum) equal to Equation 13 and solving for surface density given an assumed characteristic particle size. This surface density measurement can then be used to normalize the total disk surface density profile, thus giving us an independent estimate that can be compared to the profile derived via the disk dust lines.

While this method requires the assumption of a characteristic particle size, we note that a wide range of particle sizes will result in the same drift ice line location (see Figure 5) and thus, the derived surface density profile will not be sensitively affected by this value. For computational purposes we use a particle size of 1 m such that its stopping time exceeds unity for a wide range of disk parameter space. As is the case for all particles close to or larger than a stopping time of unity, these particles should desorb after drifting to the maximum drift ice line.

Consideration of the drift ice line in determining disk surface density has the advantage of not needing a measure of the molecular abundance in the disk. Regardless of CO abundance this location will be constant as it is dependent on disk surface density and not the CO surface density. Thus, a sensitive probe of the CO ice line (i.e. the C18O $J=3-2$ emission) should be able to detect the uptick in CO emission past this point and provide a probe of disk surface density from this measurement alone without further assumption.

5.3.2 Classical Ice Line Surface Density Determination

If the disk in question is in the classical ice line regime (see Section 5.2), the surface density can be derived by balancing Equation 11 and Equation 12. Unlike the drift dominated regime, this calculation requires an assumed molecular abundance. While these uncertainties may diminish the robustness of this test, we can make the simplifying assumption that molecular abundances are roughly constant across a single stellar type. This could be a reasonable assumption if molecular abundances are primarily shifted from ISM values via photochemical processes. If we accept this assumption as true, we can determine the disk surface density profile for disks of a stellar type for which one member has a relatively well determined molecular abundance. We note that, for this method, a factor of 2 change in the ice line location will change the derived surface density by a factor of 3.

The relative behavior of the classical ice line across many disks will also improve this test and may indeed allow the potential for deriving the CO abundance in conjunction with Test 1.

• •

Test 2 Observational Requirements: This test can be carried out on disks with accurate ice line measurements (either via C18O or N2H+). To carry out this test more precisely it is also useful to have observationally derived CO surface density profiles for these disks.

5.4 Test 3: Disk Dust and Ice Line Scalings

The third test is whether or not the dust and ice lines scale oppositely across a large sample of disks. Our model predicts that this will happen as disks with dust lines at larger radial scales should have ice lines located at shorter radii across a particular disk temperature profile and molecular abundance. This arises as larger particles become well coupled to the gas as disk surface density increases, thus causing a given particle size to have a dust line located further out in the disk. Conversely, increased surface density increases the adsorption flux onto a grain such that the freezing temperature is hotter and the ice line is moved closer to the star.

To clearly see how the dust and ice lines in a disk scale as a function of surface density alone, we vary the disk surface density for a set disk temperature profile. Using the empirical evidence that the Earth is not formed of water ice, we consider the case of a passively irradiated disk that is normalized such that the H2O ice line is outside of 1 AU as shown in Figure 12 as our fiducial temperature profile. We take this temperature profile to be roughly representative of a young sun like star. This profile follows Equation 1 with a derived $T_{0}\approx 210$ K.

This profile is derived through the use of an assumed canonical H2O abundance of 0.9$\times 10^{-4}$ n${}_{\text{H}}$ (Pontoppidan, 2006). While this assumption is sufficient for our illustrative example, we again note that disk molecular abundances may be poorly constrained. However, let us naively assume that molecular abundances vary from ISM values as a function of photochemistry such that they are constant across a particular stellar type and thus a particular temperature profile. Thus, using this temperature profile we can derive the trends that we expect to see as we vary the disk surface density.

We find that for all three of our molecular species: CO, CO2, and H2O, the radial extent of the disk increases with increasing surface density while the ice line location decreases as expected (see Figure 13). For the CO and H2O ice lines, the drift distance past the classical ice line is negligible while for the CO2 ice line drift plays a small role that decreases in importance with increasing surface density. We note that while the importance of the drift ice line will vary with disk parameters, the general trend should hold.

This trend remains true across a range of molecular abundances as long as they are held constant across a sample of disks. Here we have assumed the following abundance values for our molecular species: $\text{n}_{\text{CO}}=1.5\times 10^{-4}\;\text{n}_{\text{H}}$, $\text{n}_{\text{CO}_{2}}=0.3\times 10^{-4}\;\text{n}_{\text{H}}$, $\text{n}_{\text{H}_{2}\text{O}}=0.9\times 10^{-4}\;\text{n}_{\text{H}}$ (Pontoppidan, 2006). We note that lower molecular abundances result in lower freezing temperatures (see Section 4.1) which, depending on the temperature profile of the disk, pushes the classical ice lines further from the star with the reverse being true for higher molecular abundances.

This trend is therefore a key observational diagnostic that will be confirmed if, across a constant molecular abundance, disks with dust lines at larger radii have ice lines closer to their star. We also note that, if drift determines the ice line location (see Section 5.2) we can observe these trends without the need for an assumed molecular abundance.

• •

Test 3 Observational Requirements: To best perform this test there needs to exist a significant sample of disks with observations of both the dust and ice lines.

6 Summary & Discussion

We propose a novel method to derive disk surface density through the consideration of disk dust and ice lines. To derive this method we adopt an agnostic point of view in regards to disk surface density which we do in particular response to the uncertainties that accompany typical observational tracers. This method relies on the assumption that, at late stages of evolution, the growth timescale, drift timescale, and the lifetime of the disk are all equal for the dominant particles at a dust line location ($s=\lambda_{\text{obs}}$).

While other work finds that these timescales are equal they are often unable to match the disk surface density. We therefore make the assumption that these timescales are equal and use this to determine the surface density profile without evoking other observational tracers. These assumptions allow us to self-consistently derive a disk surface density profile as well as a dust-to-gas ratio in the outer disk without the further assumption of a given molecular abundance.

We apply our modeling technique to our fiducial disk TW Hya. We find that our derived surface density profile and dust-to-gas ratio are consistent with the lower limits found through measurements of the HD gas in the disk (Bergin et al., 2013; Schwarz et al., 2016). Using our derived surface density profile we uncover a theoretical estimate of the disk accretion rate that is more closely aligned with the measured accretion onto TW Hya. We further find that our theoretical classical and drift ice lines have clear observational analogues where the classical ice line is predictive of the ice line found via N2H+ emission (Qi et al., 2013) and the drift ice line is probed by the C18O emission (Nomura et al., 2016). We conclude that the ice line derived through observations of C18O emission may be more sensitive to the extent of the drift ice line in the disk and thus may be the best observational method to test our model’s assumptions. We note here that our method highlights the likelihood of detecting multiple ice lines for a given molecular species as there will be relatively large regions in the disk where some of the material has desorbed and some of the material has not.

Furthermore, if the banded structure in the ALMA images of HL Tau and TW Hydra reflect substantial variations in surface density then our model should still hold. However, ice lines may be preferentially found in less dense bands in the disk where drift slows.

We next consider a large range of theoretical disk parameter space and uncover three observational tests of our model. The first test is whether or not the dust lines of other disks, once converted to surface densities, can be matched with a previously derived normalized surface density profile for the disk. For TW Hya, the disk dust lines fall near the exponential cut off region of the disk but there is no reason that this needs to be true across many objects and thus provides a powerful test of the model. In the age of ALMA, dust lines can be determined observationally for a significant sample to protoplanetary disks which will provide the ideal observations necessary for this test.

The second test is whether we can derive a surface density profile from the disk ice lines that matches the surface density profile derived via disk dust lines. The third and final test is whether or not disk dust and ice lines scale oppositely, as a function of surface density, across a large sample of disks. Disk ice lines have been observed for an increasing number of disks such that the last two tests could be carried out in the near future with the aid of facilities such as ALMA.

7 Acknowledgements

This work was done as a part of the Kavli Summer Program in Astrophysics: Exoplanetary Atmospheres 2016. We would like to thank Jonathan Fortney and Pascale Garaud for providing us with the opportunity to work together on this project and for organizing a wonderful workshop. We thank James Owen, Michael Rosenthal, Catherine Espaillat and Laura Pérez for their useful comments and insightful discussion. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant DGE1339067. RMC acknowledges support from NSF grant number AST-1555385.