Disk masses around solar-mass stars are underestimated by CO observations

TL;DR

This paper demonstrates that CO observations significantly underestimate the true gas mass in protostellar disks due to chemical depletion and optical thickness, suggesting improved observational strategies for accurate mass measurements.

Contribution

The study introduces radiative transfer simulations showing how CO depletion and optical thickness affect disk mass estimates, proposing new observational strategies for accuracy.

Findings

CO-based mass estimates are underestimated by over an order of magnitude.

CO abundance varies with disk radius and time, causing age/mass degeneracy.

Recommended strategies include using low-J transitions and multiple isotopologues.

Abstract

Gas in protostellar disks provides the raw material for giant planet formation and controls the dynamics of the planetesimal-building dust grains. Accurate gas mass measurements help map the observed properties of planet-forming disks onto the formation environments of known exoplanets. Rare isotopologues of carbon monoxide (CO) have been used as gas mass tracers for disks in the Lupus star-forming region, with an assumed interstellar CO/H$_2$ abundance ratio. Unfortunately, observations of T-Tauri disks show that CO abundance is not interstellar---a finding reproduced by models that show CO abundance decreasing both with distance from the star and as a function of time. Here we present radiative transfer simulations that assess the accuracy of CO-based disk mass measurements. We find that the combination of CO chemical depletion in the outer disk and optically thick emission from the inner disk leads observers to underestimate gas mass by more than an order of magnitude if they use the standard assumptions of interstellar CO/H$_2$ ratio and optically thin emission. Furthermore, CO abundance changes on million-year timescales, introducing an age/mass degeneracy into observations. To reach factor of a few accuracy for CO-based disk mass measurements, we suggest that observers and modelers adopt the following strategies: (1) select the low-$J$ transitions; (2) observe multiple CO isotopologues and use either intensity ratios or normalized line profiles to diagnose CO chemical depletion; and (3) use spatially resolved observations to measure the CO abundance distribution.

Equations

$$\mathcal{N}=\frac{\mathcal{N}_{J}}{g_{J}}\times Q\times e^{\mbox{$E_{\rm up}$}/T}$$

$$M=\frac{\mathcal{N}\mbox{$\mu_{\rm H}$}\mbox{$m_{\rm H}$}\mbox{$f_{\rm iso}$}}{\mbox{$f_{\rm C}$}\mbox{$f_{\rm CO}$}}=1.668\times 10^{-53}\mbox{M${}_{\odot}$}\ \mathcal{N}\mbox{$f_{\rm iso}$}/\mbox{$f_{\rm CO}$}$$

$$\frac{{\mathcal{N}}_J}{g_J}=\frac{L_J}{g_J{A}_{J}h\nu }$$

$$L_{J}=F_{l}({\rm cgs})\times 4\pi D_{\rm cm}^{2}=1.1964\times 10^{20}\times\mbox{$F_{l}({\rm Jy-km\ s^{-1}})$}\ (\nu/c)D_{\rm pc}^{2},$$

$$A_{J}=\frac{64\pi^{4}}{3h}\times(\nu/c)^{3}\times|\mu(J,J-1)|^{2}=3.13613\times 10^{-7}\times(\nu/c)^{3}\frac{J}{(2J+1)}\mbox{$\mu_{\rm D}$}^{2}$$

$$\frac{\mathcal{N}_{J}}{g_{J}}=6.468\times 10^{45}\frac{\mbox{$F_{l}({\rm Jy-km\ s^{-1}})$}\mbox{$D_{\rm pc}$}^{2}}{\mbox{$B_{\rm GHz}$}^{3}J^{4}\mbox{$\mu_{\rm D}$}^{2}}$$

Full text

††software: RADMC (http://www.ita.uni-heidelberg.de/ dullemond/software/radmc-3d/), LIME (Brinch & Hogerheijde, 2010)22affiliationtext: Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon, 34055, Korea

1 Introduction

Solar System formation models that allow the giant planets to form within observed protostellar disk lifetimes of a few million years (Haisch et al., 2001) often require density enhancements up to an order of magnitude above the minimum-mass solar nebula (MMSN) (Pollack et al., 1996; Hubickyj et al., 2005; Thommes et al., 2008; Lissauer et al., 2009; Dodson-Robinson & Bodenheimer, 2010; D’Angelo et al., 2014). Indeed, planet accretion may be a fundamentally inefficient process, with both collisional fragmentation (e.g., Stewart & Leinhardt, 2012) and planetesimal scattering (e.g., Ida & Lin, 2004) contributing to mass loss during solid embryo growth. Yet disk masses inferred from dust emission in (sub)millimeter often do not reach the MMSN mass of $0.01M_{\odot}$ (Weidenschilling, 1977; Hayashi, 1981), and are more commonly of order 1-10 Jupiter masses (Andrews & Williams, 2007; Williams & Cieza, 2011). Reporting on a survey of T-Tauri stars in Lupus, Ansdell et al. (2016) suggested that 80% of the disks had dust-derived total masses of $<0.01M_{\odot}$. However, dust-based disk mass estimates may be systematically low, because dust continuum observations lose sensitivity to solids that are much larger than the observing wavelength (Williams & Cieza, 2011). In addition, the standard assumptions that the dust has a single temperature and that the sub-mm emission is optically thin everywhere in the disk may not be correct. Finally, gas masses may not be related to the dust masses by the usual interstellar ratio of 100. It is essential to measure the gas mass of disks directly.

One such gas mass measurement came from Bergin et al. (2013), who used Herschel observations of the HD $(J=1\rightarrow 0)$ transition to calculate a mass of $0.06M_{\odot}$, or 6 MMSN, for the disk surrounding TW Hydra—a surprisingly large mass given the star age of $\sim 10$ Myr. The HD lines have since been detected in two more disks (McClure et al., 2016), and are consistent with gas masses of 1-4.7 MMSN (DM Tau) and 2.5-20.4 MMSN (GM Aur). Although the GREAT instrument on the far-IR observatory SOFIA111https://www.sofia.usra.edu/ covers the frequency of the HD transition, it is not sensitive enough to observe HD in nearby disks. In the absence of the capability to observe the HD $J=1\rightarrow 0$ transition in disks, CO has been the standard tracer of the gas mass because it is believed to have simple chemistry and to stay in the gas phase wherever $T>20$ K in disks around Sunlike stars (Öberg et al., 2011; Qi et al., 2013), a region that includes the entire vertical column in the inner 30 AU and the warm surface layers of the outer disk. While the emission from 12C16O is typically optically thick, it has been suggested (van Zadelhoff et al., 2001; Dartois et al., 2003) that rare isotopologues of CO could be used to probe the disk midplanes. Yu et al. (2016) (hereafter Paper 1) have shown that the vertical optical depth of low-J rotational emission lines of C17O is around unity in the inner $\sim 20$ AU of a 1.5-MMSN disk, meaning observers could see emission from the disk midplane—where most of the mass is concentrated—using C17O lines. Unfortunately, Paper 1 also revealed some complexities in the CO chemistry that would interfere with disk mass measurements. First, the CO abundance varies with distance from the star within the planet-forming region. Second, the CO/H2 ratio drops to an order of magnitude below the interstellar value well inside the CO freeze-out radius because of chemical depletion of CO. Finally, the CO abundance is a function of time, which introduces an age-mass degeneracy into the interpretation of the observations.

Recent attempts to calculate gas-to-dust mass ratios using observations of rare CO isotopologues have also revealed problems. In their survey of disks in the Lupus star-forming region, Ansdell et al. (2016) found gas-to-dust ratios, calculated assuming a constant CO/H2 ratio of 10*-4*, to be much lower than the interstellar value of 100. Yet the stars in the Lupus sample are still accreting, indicating that abundant gas is present. Studies of a more massive disk and star (2.3 M⊙) also found very low gas-to-dust ratios (Boneberg et al., 2016). Correcting for the joint effects of freeze-out and isotope-selective photodissociation still does not bring the estimated gas-to-dust ratios up to 100, according to the chemical models of Miotello et al. (2017). Either the disks are in the process of dispersing—unlikely given the rapid depletion timescale of $\sim 10^{5}$ years once photoevaporation dominates the disk dynamics (e.g., Hollenbach et al., 2000; Alexander et al., 2014; Gorti et al., 2016)—or additional pathways to remove gaseous CO (hereafter referred to as chemical depletion) become important in disks, as suggested by Dutrey et al. (2003), Favre et al. (2013), Miotello et al. (2017), and McClure et al. (2016). Some disks around HAeBe stars also appear to have CO/H${}_{2}<10^{-4}$ (e.g. Chapillon et al., 2008; Bruderer et al., 2012).

Also challenging the assumption of simple CO chemistry and a constant CO/H2 ratio in regions where CO is not frozen out are the TW Hya observations of Schwarz et al. (2016), which show a drop in CO column density at $\sim 20$ AU, while the dust column density remains roughly constant with radius. In Paper 1, we found that CO chemical depletion due to dissociation by He+ and subsequent complex organic molecule (COM) formation causes the CO/H2 abundance ratio to drop far below the interstellar value of $10^{-4}$ at $r>20\,{\rm AU}$, well inside the CO freeze-out radius in the model (see figure 1). For disks with similar chemistry to the Paper 1 models, using interstellar abundance ratios to extrapolate from CO to H2 column density would result in large underestimates of disk mass.

In this work we use radiative transfer models of CO rotational emission to assess the usefulness of rare CO isotopologues as disk mass indicators, given the possibility of chemical CO chemical depletion. We summarize the key results from our chemical evolution models in section 2, including a new model for a more massive disk. We present the setup for and results of our line radiative transfer models in section 3. In section 4, we demonstrate that standard gas-mass measurement methods—based on integrated fluxes of lines assumed to be optically thin—fail when applied to the simulated CO emission from our model disk. We also evaluate the performance of published optical-depth correction methods (Williams & Best, 2014; Miotello et al., 2016), showing that they are inadequate for disks with non-uniform CO abundance. In section 5 we further highlight the age-mass degeneracy problem caused by CO chemical depletion. Next we show that (1) combining observations of multiple isotopologues, (2) using information on line profiles, and (3) examining the spatial distribution of CO can diagnose CO chemical depletion (section 6). In section 7, we explore whether our predicted CO/H2 abundance ratios can increase observed masses up to MMSN or higher when applied to the Lupus disk sample of Ansdell et al. (2016). Section 8 summarizes our results.

2 Disk Model

We adopt the chemical-dynamical model from Paper 1 as the basis for our line radiative transfer models. Paper 1 presented the chemical evolution of a $0.015\,{M_{\odot}}$ disk around a Solar-type star for $3$ Myr. In this paper, we add a chemical evolution model for a disk that is twice as massive, at $0.03\,{M_{\odot}}$. The mass distributions and accretion temperatures of both model disks were presented by Landry et al. (2013), and we find the stellar contribution to the disk heating from the dust radiative transfer code RADMC222http://www.ita.uni-heidelberg.de/ dullemond/software/radmc-3d/; developed by C. Dullemond. Because our goal is to measure disk masses in giant planet-forming regions, we focus our study on the region covered in Landry et al. (2013) – the inner $70$ AU of the disk. (As the disk ages, the amount of mass in the inner $70$ AU decreases due to viscous spearing and accretion – all disk masses reported in this paper are calculated as $\int 2\pi r\Sigma~{}dr$, where $\Sigma$ is the column density, out to 70 AU.) The chemical reaction network is run locally at each independent (r, z) grid point for $3$ Myr, under the assumption that the chemical reaction timescale is much shorter than the viscous timescale—an assumption that is true for freezeout, desorption, and grain-surface reactions, but which may fail for gas-phase reactions. Each disk gridpoint starts with gas and ice abundances resulting from a 1 Myr simulation of the chemical evolution of a parent molecular cloud; as a result, a substantial fraction of the carbon is tied up in CO2 and other ices at the start of disk evolution.

The chemical evolution models include C, H, O, N based on the UMIST database RATE06 (Woodall et al., 2007). Woods & Willacy (2009) extended the network to include C isotopes, and we included both C and O isotopes in Paper I. The chemical models follow the chemistry of $588$ species, $414$ gas-phase and $174$ ices, for $3\,{\rm Myr}$ from the beginning of the T-Tauri phase. The reaction network contains gas-phase reactions (including those that lead to C and O fractionation), grain-surface reactions, freezeout, thermal desorption, and reactions triggered by UV, X-rays and cosmic rays, such as isotope-selective photodissociation. To make our simulations computationally tractable, we include only molecules with two or fewer carbon atoms, which limits the network to 11316 reactions. Some grain-surface hydrogenation reactions, such as $(\mathrm{C}_{2}\mathrm{H}_{5}\;{\rm ice})+H\rightarrow(\mathrm{C}_{2}\mathrm{H}_{6}\;{\rm ice})$, are not included because they only lead to sinks; we can save computational time by not computing the associated reaction rate and letting $\mathrm{C}_{2}\mathrm{H}_{5}$ ice be the sink instead of $\mathrm{C}_{2}\mathrm{H}_{6}$ ice.

The simplifications we have made to the reaction network do not artificially remove carbon from the gas phase. In Paper 1 we demonstrated that using a range of chemical networks, input chemical abundances and dust properties does not change the evolution of CO abundance significantly. We take the temperatures and abundances from the fiducial model in Paper 1 as the foundation for this study and refer readers to Paper 1 for detailed discussions of disk model assumptions. The most debatable assumption from Paper 1 is that abundances in gas and icy components are inherited from the molecular cloud without modification as they enter the disk. Indeed, Drozdovskaya et al. (2016) suggest that disk midplane composition is largely determined by the conditions during cloud collapse. Visser et al. (2009) find that CO ice formed in molecular clouds desorbs during cloud collapse, though it re-freezes without significant abundance modification where the disk temperature is $<18$ K.

The disk is primarily heated passively from stellar irradiation; accretion contributes very little to the energy budget (Landry et al. 2013, Paper 1). For the temperature calculation, we use the dust opacities computed by Semenov et al. (2003), adjusted for our gas/dust ratio. We assume the dust has already grown and aggregated at the start of the T-Tauri phase of disk evolution (Oliveira et al., 2010; Pérez et al., 2012; Birnstiel et al., 2011; Garaud et al., 2013), and 90% of the solids have grown to still larger sizes (pebbles, rocks, etc.), that contribute very little submillmeter emission. We distinguish these solids from dust and take a gas/dust mass ratio of 1000. In Paper 1, we demonstrated that the CO abundance is only weakly sensitive to evolution of the gas/dust ratio from 100 to 1000. Our grain model assumes that there is a constant replenishment of small grains by collisions between larger objects (see Dullemond & Dominik (2005); Brauer et al. (2008); Birnstiel et al. (2009); Wada et al. (2008, 2009); Windmark et al. (2012); though Zsom et al. (2010) argues against replenishing micron-size dust grains through collisions) and the size distribution of dust does not evolve. The theory that most of the solid mass has aggregated into pebbles and larger objects invisible to submillimeter wave observations is the second major, non-standard assumption in our model (though Andrews et al. (2012) and Cleeves et al. (2015) modeled the TW Hya disk with two distinct grain populations, small grains that provide the visible/infrared opacity and large grains with a maximum size of 1 mm that provide the submillimeter opacity). Following Cleeves et al. (2015), we assume the 90% of solid mass that has grown to pebble and larger sizes provides negligible surface area for chemical reactions and exclude it from the chemical model.

As the central star dims during its time on the Hayashi track, the disk cools and its scale height shrinks. The computational surface of our model grid, defined as the layer where the optical depth to the disk’s own radiation is $\tau=0.2$ (Landry et al., 2013), moves from about two pressure scale heights to about one scale height above the midplane as the disk becomes cooler and thinner over time. The chemical model therefore contains fewer grid cells at the end of evolution than at the beginning, as grid layers high above the midplane begin to empty out. In appendix A, we show that our line profile and intensity calculations are only weakly sensitive to the changing disk surface.

Stellar heating is efficient enough to prevent CO from freezing out in our modeled region—the inner 70 AU of the disk—at all vertical heights and at any time in the $3$ Myr of evolution. However, CO is depleted beyond $20$ AU from the central star due to the formation of complex organic molecules. The CO chemical depletion is driven by ionization of helium from X-rays and cosmic rays and happens over a million-year time scale. As a result, the CO abundance changes both with location in the disk and with time. We show the color map of CO abundance at $2$ Myr of the disk evolution in Fig. 1 for both the Paper 1 disk ($0.015$ M⊙; left) and the new model of the more massive disk ($0.03$ M⊙; right). The abundance is defined as the ratio of the number density of CO with respect to the number density of hydrogen nuclei (n${}_{\rm H}+2\rm n_{\rm H_{2}}$). After 2 Myr, the CO chemical depletion front moves inward to 20 AU, with the surface layers more depleted than the disk midplane. Fig. 2 shows the abundances of major carbon-bearing species as a function of time for the midplane at 40 AU in each disk, demonstrating the gradual chemical depletion of CO and subsequent sequestration of carbon into organic ices. We note that our network is not extensive enough to determine the exact end-product of organic ice formation. The most complex hydrocarbon we include is C2H5, which in reality should hydrogenate to form ethane, and radicals of the form C2Hx should react with carbon atoms or hydrides to form longer carbon chains. However, these complex products of organic grain surface chemistry are all less volatile than the C2Hx species and would stay on the grain surfaces unless the grains drift radially inward (e.g. Birnstiel & Andrews 2014) or experience a transient heating event (e.g. Cody et al. 2017, Cieza et al. 2016, Vorobyov & Basu 2015). Since almost all resolved T-Tauri disks detected in dust continuum emission have radii much larger than the $\sim 40$ AU C2H5 ice line (Paper 1) in our models333See the catalog of resolved disk images at www.circumstellardisks.org, we do not expect radial drift to deposit a significant amount of organic gas that could be recycled to form CO in the inner disk. To the extent that there are real astrophysical disks that evolve quiescently during the T-Tauri phase, our conclusion that carbon liberated by CO chemical depletion becomes sequestered in ices is robust. The more massive disk shows a similar pattern of CO chemical depletion to the disk presented in Paper 1, though the depletion timescale is a bit longer because the higher column density decreases the ionization fraction—and thus the abundance of ionized helium—in the midplane.

The net result of all chemical models presented here and in Paper 1 is that CO becomes severely depleted well inside the CO freeze-out radius in disks with masses above the minimum needed to form planetary systems. Similar effects have been seen in other chemical evolution calculations (Aikawa et al. 1997, 1999; Furuya & Aikawa 2014; Walsh et al. 2014; Bergin et al. 2014). Aikawa et al. (1997) and Aikawa et al. (1999) first pointed out that CO can react to form less volatile molecules such as CO2, HCN, H2CO, CH4 and larger hydrocarbons over Myr timescales. Aikawa et al. (1997) used a static disk (not evolving), and Aikawa et al. (1999) used a vertically isothermal model and found a larger number of simple molecules as products (as opposed to a few complex molecules), but the essence of the chemical depletion of CO is the same as found here and in Paper 1. Walsh et al. (2014) simulated the composition of complex organic molecules in a disk with no temperature evolution for about $1$ Myr and found the formation of complex organic molecules in the disk midplane via grain-surface reactions, while Bergin et al. (2014) found gas-phase organic formation and subsequent freezeout onto grain surfaces. Furuya & Aikawa (2014) investigated the carbon and nitrogen chemistry during turbulent mixing, and found that volatile transport enhances COM formation near the surface and suppresses it in the disk midplane. It is difficult to compare our results directly with the earlier calculations because our models have different central star properties and less accretion heating, but it is clear that a variety of chemical processes can force a disk’s CO/H2 abundance ratio far from the interstellar value.

In Figure 3 we show the fraction of our model disk’s carbon atoms contained in CO gas ($f_{\rm CO}$) as a function of time. The fraction is initially low because evolution in the molecular cloud has sequestered carbon in ices; for a while these evaporate, increasing $f_{\rm CO}$, but then CO chemical depletion and formation of icy organics cause $f_{\rm CO}$ to decrease. To complicate matters, the abundance of CO actually increases with time at small radii, as CO2 ice desorbs and dissociates to form CO gas. Although the $0.015\,{M_{\odot}}$ disk suffers more CO chemical depletion at larger radii, the re-formation of CO gas in the inner disk also proceeds at a higher rate than in the $0.03\,{M_{\odot}}$ disk, leading to a higher disk-averaged $f_{\rm CO}$ for the $0.015\,{M_{\odot}}$ disk. These facts will compromise attempts to measure disk masses using CO isotopes (§4). Worse yet, the CO abundance is a function of time, leading to an age-mass degeneracy in interpreting observations (§5).

3 Line radiative transfer model and mass estimates

In order to understand the effects of chemical evolution on CO emission, we set up radiative transfer models with the publicly available code LIME (LIne Modeling Engine; Brinch & Hogerheijde, 2010). LIME calculates either non-LTE or LTE line excitation and solves the radiative transfer problem for molecular gas in arbitrary 3D geometries. For recent examples of studies using LIME, see Walsh et al. (2016) and Öberg et al. (2015). We adopt the energy levels from the Leiden Atomic and Molecular Database (LAMDA)444http://home.strw.leidenuniv.nl/ moldata/. As a first order approximation, we do not consider the hyperfine splitting in C17O emission.

As in Paper 1, we model emission from within a radius of $70\,{\rm AU}$ of the central star, which corresponds to a $1\arcsec$ beam diameter for an assumed distance of $140\,{\rm pc}$ from the Sun. We consider line broadening due to Keplerian rotation, thermal velocity and micro-turbulence. Thermal velocities are calculated assuming a Maxwell-Boltzmann speed distribution based on the disk’s temperature structure from Paper 1. For the micro-turbulence, we again assume an isotropic Maxwell-Boltzmann speed distribution with RMS of $100\,{\rm m/s}$, consistent with the upper limit to microturbulent speed found by Flaherty et al. (2015) in a fit to multiple CO emission lines in the HD 163296 disk. For computing the level populations, we set a minimum scale of $0.07\,{\rm AU}$ to guarantee sub-pixel sampling of both Keplerian speeds and CO abundance gradients. We first generate the synthetic datacube of intensity as a function of $x$, $y$, and velocity for a disk around a $0.95\,{M_{\odot}}$ star at $30\arcdeg$ inclination, similar to the disk surrounding AS 209. In velocity space, the spectra have $300$ channels of $125\,{\rm m/s}$ resolution. At any specific velocity, the synthetic image contains $600\times 600$ pixels of $0.003\arcsec\times 0.003\arcsec$ in size. Finally, we generate the synthetic spectra presented here by integrating each velocity component over a square with $1.2\arcsec$ sides ($400\times 400$ pixels), larger than the angular size of the disk. The pixels not covered by the disk contribute no flux and are included simply for ease of integration—here we assume that the sky background contains negligible flux compared with the disk at all wavelengths.

Our current models assume that the gas temperature is the same as the dust temperature. This is a valid assumption in estimating the optical depth of C17O and C18O, as done in Paper 1, because the emission primarily comes from disk midplane and the midplane is strongly shielded from UV radiation. However, hot gas on the disk surface is more emissive, and we would need to consider the difference between the gas and dust temperature in order to use our models to fit high-J spectral lines emitted from the disk surface. Similarly, we would need to revisit the temperature calculation for a disk surrounding a star with a stronger UV field, which could decouple the gas and dust temperatures. In Appendix A, we demonstrate that the LTE approximation is adequate for computing the energy level populations, and show that the decreasing disk scale height has little effect on the computed line profiles (see §2). For our purposes, we have adequately modeled the CO rare isotopologue emission, even though our models do not extend vertically to a large number of scale heights. For the rest of the paper, we focus our line profile discussion on $J=3\rightarrow 2$ and $J=2\rightarrow 1$, which are the most commonly observed transitions.

Figure 4 shows one set of results from our line radiative transfer models—the time evolution of the $J=3\rightarrow 2$ line profiles for three isotopologues and two model disks. Line profiles in Janskys are plotted as well as normalized profiles, which better show the evolving shape of the lines (see §6.2 for more on diagnosing CO chemical depletion by comparing line profiles from multiple isotopologues). Figure 5 shows the total intensities (integrated in velocity space) of $J=3\rightarrow 2$ emission from all CO isotopologues as a function of time. The decline with time is quite dramatic, especially for the rarer isotopologues where it exceeds an order of magnitude over 3 Myr. The emission from rarer isotopologues is very weak at later times, which can be expensive to observe, especially if one wants to achieve enough signal-to-noise to study the line profiles.

Finally, we use standard observational procedures to try to recover the apparent masses of the disks, as an observer would do. Here we will use the $J=2\rightarrow 1$ line, as it is commonly used to obtain disk masses (e.g.,Williams & Best 2014, though Ansdell et al. (2016) use the $J=3\rightarrow 2$ line). We show in Appendix C that higher $J$ transitions underestimate the mass more severely than $J=2\rightarrow 1$ and $J=3\rightarrow 2$. From our model line profiles we calculate level populations, total number of CO molecules, and finally gas mass. The equations used to derive the number of CO molecules in each energy level ($\mathcal{N}_{J}$) from integrated line flux, assuming optically thin emission, are given in Appendix B. Values of $\mathcal{N}_{J}$ can then be translated to mass with the following equations. The total number of CO molecules in the disk $\mathcal{N}$ is given by

[TABLE]

where $E_{\rm up}$ is the upper-state energy of the transition in K and $Q$ is the partition function, generally approximated by $Q=kT/hB$ (with $B$ in Hz). Note that the temperature of the CO reservoir appears in both the exponential and the partition function in Equation 1; in Appendix C.1 we show that no single temperature describes the CO reservoir. We then compute the mass of gas from

[TABLE]

where $\mu_{\rm H}$ is the mean atomic weight of the ISM including He, $m_{\rm H}$ is the mass of a hydrogen atom, $f_{\rm iso}$ is the ratio of the mass of the most common isotopologue (12C16O) to that of the one being used, $f_{\rm C}$ is the abundance ratio of carbon to hydrogen nuclei, and $f_{\rm CO}$ is the fraction of C in CO, averaged over the disk. \mbox{\mu_{\rm H}}=1.43 and our chemical model has \mbox{f_{\rm C}}=7.21\times 10^{-5}.

In the next section we describe how equations 1 and 2 to obtain the correct mass of our model disk.

4 Even rare CO isotopologues underestimate disk mass

12C16O and 13C16O will be optically thick in disks, so the usual approach is to observe rarer isotopologues. At first glance, the CO chemical depletion we predict might be expected to ameliorate optical depth problems, but the concentration of CO to small radii increases optical depths there. At the same time, the chemical depletion of CO at larger radii causes mass underestimates by lowering the value of $f_{\rm CO}$. The net result is that simple analysis fails for the following reasons.

First, the usual assumption is that \mbox{f_{\rm CO}}=1 inside the CO ice line, but our models challenge that assumption. Figure 3 shows how $f_{\rm CO}$ varies with time, never rising above 0.3 ($0.015$ M⊙ disk) or $0.25$ ($0.03$ M⊙ disk) and approaching 0.12 by 3 Myr. Since our disk does not include radii where CO itself would freeze, the low $f_{\rm CO}$ is due to what we call chemical depletion. Alone, with no other error sources in equations 1 and 2, it will cause underestimates of disk mass by factors of 3 to 8.

Second, the concentration of CO toward the inner disk means that temperatures, and hence partition functions, are higher than usually assumed for much of the CO reservoir. The most common assumption for temperature is $T=20$ K (e.g., Ansdell et al. 2016). A temperature averaged over the the model disk and weighted by the number density of CO molecules ranges from 50 to 70 K, decreasing as the star and disk evolve, but stabilizing around 55 K for the 0.015 M⊙ disk after about 1 Myr because the increasing concentration of CO in the inner disk counteracts the dropping luminosity of the star (Figure 6). If CO chemical depletion is operating in the outer disk, assuming $T=20$ K for the gas where CO is concentrated is never a good choice; it will systematically underestimate masses by a factor of about 2. In disks around T-Tauri stars that are less luminous than the proto-sun at 3 Myr, or that receive a low cosmic-ray flux so that cosmic ray-induced photons do not generate CO from CO2 gas in the inner disk, 20 K may be appropriate, but we recommend constructing model-based temperature estimates for the CO reservoir rather than making any assumptions. As a further complication, we show in Appendix C that rotation diagrams from multiple transitions are not effective in determining the best temperature to assume.

Third, the increase in CO abundance in the inner disk causes optical depth effects where most of the CO actually resides (see Fig. 10 of Paper 1). Even if we correct for the mean $f_{\rm CO}$ and somehow find an accurate value of $T$ averaged over the CO reservoir, the standard analysis still underestimates the mass. Figure 7 shows the mass estimates from different isotopes and analysis procedures versus time, along with the actual disk mass, which decreases slightly due to disk spreading and accretion onto the star. The estimates with black symbols are computed with equations 1 and 2 and based on the assumption of optically thin emission in the $J=2\rightarrow 1$ line of the listed isotopologues. They all use the correct $f_{\rm CO}$ and the mass-weighted CO temperature for our model—information that would not be available for astronomical sources—yet they all still underestimate the gas mass substantially. The rarer isotopologues perform best, suggesting that optical depth is the main culprit. The recent analysis of the Lupus disks (Ansdell et al., 2016) by Miotello et al. (2017) also suggests optical depth effects, even after ratios of 13CO and C18O were used to correct for optical depth. The apparent gas to dust ratio declines with increasing dust mass in their Figure 6; if the dust mass is a reasonable proxy for disk mass, a more massive disk will have larger (and more complicated) effects from optical depth. Our results indicate that models that do not account for the radial dependence of $f_{\rm CO}$ will underestimate the gas mass more severely for more massive disks, as observed in Figure 6 of Miotello et al. (2017).

4.1 Correcting for optical depth using multiple isotopologues

One approach to dealing with optical depth is to observe the same rotational transition in multiple isotopologues (e.g., Williams & Best, 2014; Miotello et al., 2016). We tried a simple correction for our 13CO-based mass estimate by comparing with C18O, using the C18O/13CO line intensity ratios plus the initial isotope ratios from the chemical model. This correction yielded the hollow green hexagons in Figure 7, which still underestimate disk mass badly, performing only as well as C18O observations would on their own.

More sophisticated optical depth corrections can be made with fitting formulae that relate the $J=2\rightarrow 1$ line intensity of either 13CO or C18O (or their ratio) to the total disk mass, based on a suite of models (e.g., Williams & Best, 2014; Miotello et al., 2016). We placed our model line luminosities on Figure 7 (left panel) of the Williams and Best models to estimate mass, producing the red points, which have no correction for $f_{\rm CO}$. They do well for early times when the CO is more uniformly distributed over the disk, but underestimate the mass badly at later times when $f_{\rm CO}$ drops especially at larger radii. Most of our points lie near the top of their distribution of models, and interpolation between mass models is uncertain by a factor of three. Second, we used the formulae in equation 2 with coefficients in Table A.1 of Miotello et al. (2016) to estimate disk mass from luminosities of the $J=2\rightarrow 1$ lines of 13CO and C18O. Those from 13CO underestimated the mass badly, but those from C18O did better (blue octagons). In these comparisons, we adjusted their mass estimates to be consistent with our assumption about the atomic carbon abundance, but not for the fraction of carbon in CO. In a recent paper, Miotello et al. (2017) find the same result; applying their models to the Lupus data (Ansdell et al., 2016) data leads to very low gas to dust ratios (less than 10 in most cases) unless CO is depleted. The Miotello et al. grid of chemical models was optimized to treat CO self-shielding and isotopic fractionation, which our chemical models also include; the low gas/dust ratios implied by C18O and 13CO observations must result from other CO-depletion pathways.

Clearly, both optical depth and CO chemical depletion must be accounted for when measuring disk masses: Figure 7 shows that measurements incorporating one correction, but not the other, will fail. For masses based on a single emission line, so that no optical depth correction is possible, the isotopotologue that performs the best is C17O, which underestimates the mass by a factor of 2-3 if the correct value of $f_{\rm CO}$ is known (which it will not be in real astrophysical situations). In Appendix C, we show that the $J=1\rightarrow 0$ line does much better than the more commonly used $J=2\rightarrow 1$ or $J=3\rightarrow 2$ lines. We believe that this is primarily due to its lower optical depth and lesser sensitivity to temperature.

The combination of the models of Miotello et al. (2016) and our calculations of $f_{\rm CO}$ would improve the mass estimates. For example, if we use our “insider information” (i.e. model-derived knowledge) about $f_{\rm CO}$ as a function of time for our 0.015 M⊙ disk to further correct the masses from the formulae in Miotello et al. (2016), the resulting masses are accurate to within a factor of two at all ages. Unfortunately, to apply this method to observations, one needs to know the disk age. We explore the issues arising from that fact in the next section.

5 The Age-Mass Degeneracy

The time dependence of the fraction of C in CO (see Fig. 3) implies that one needs to know the disk age to apply a mass correction factor. Ages of host stars are not generally known to better than $\pm 1$ Myr. Worse yet, the “age” in our models is in some sense a chemical age because the speed of chemical evolution will depend on the ionization sources. EXor outbursts and X-ray flares may also introduce brief periods of intense ionization (e.g. Audard et al., 2014). If we don’t know the chemical age, there is an age-mass degeneracy—a massive, old disk may look similar to a young, less massive one in CO rare isotopologue emission. Figure 8 shows an example of the age-mass degeneracy. The more massive disk reaches the same C18O/13CO total intensity ratio about 0.5 Myr later than the less massive disk for both $J=3\rightarrow 2$ and $J=2\rightarrow 1$.

CO-based mass measurements are further complicated by the fact that the depletion occurs in the outer disk, where most of the mass resides. The decline in CO isotopologue line intensity as a function of time is precipitous: in the $0.015\,{M_{\odot}}$ model, the intensity of 13CO J=3-2 emission drops from 2.72\mbox{{\rm Jy-km\ s^{-1}}} by $82\%$ to 0.48\mbox{{\rm Jy-km\ s^{-1}}} over the $3$ Myr disk evolution, and the intensity of C17O J=3-2 drops by $95\%$ from 1.15\mbox{{\rm Jy-km\ s^{-1}}} to 0.06\mbox{{\rm Jy-km\ s^{-1}}} (Fig. 5). For comparison, the total disk mass within 70 AU of the star drops by about 12% over 3 Myr and $f_{\rm CO}$ drops by a factor of 2-3. Clearly, a simple correction for $f_{\rm CO}$ does not capture all of the decrease in line emission with time. Without a good model of CO abundance as a function of radius and time, line-intensity measurements yield disk mass estimates that are accurate to factor of two at best and are systematically underestimated.

6 Diagnosing CO chemical depletion

In section 2 we demonstrated that the CO abundance decreases with both radius and time. In Section 4, we showed that other effects add to the underestimation of disk mass and that the underestimate gets worse with age. Given that star ages can be uncertain by over 1 Myr (e.g., Soderblom et al., 2014)—the timescale over which we observe CO depletion in our model—interpreting CO observations requires a more direct CO chemical depletion indicator than the star age. The need for a disk-based CO chemical depletion indicator is especially evident given the likelihood that disks with different masses or incident cosmic-ray fluxes may evolve at different speeds.

We explore here three types of observations that can diagnose CO chemical depletion . First, isotopologue intensity ratios for a transition change with time, revealing departures from the ISM CO/H2 ratio. Second, line profile shapes for the most optically thin isotopologues, C17O and C18O, widen over time as the outer disk loses CO, while 13CO and CO profile shapes evolve very little. Third, spatially resolved observations will reveal CO chemical depletion patterns. While using intensity ratios to diagnose CO chemical depletion requires less observing time, high signal-to-noise line profiles or spatially resolved observations contain valuable information that could allow observers to reconstruct the radial distribution of CO gas. We present all three strategies here, beginning with intensity ratios.

6.1 How intensity ratios reveal CO chemical depletion

The relative strengths of emission lines for different isotopologues change significantly over time. This is because different isotopologues react to the change in CO abundance differently due to their various optical depths. Re-examining Figure 8, we see that after $1\,{\rm Myr}$, the C18O/13CO intensity ratios decrease steadily as the disk-averaged CO abundance declines monotonically (Fig. 3). One could use Figure 8 to estimate the chemical age of an observed disk. However, the predicted intensity ratios are higher for the more massive disk, so there is still a 0.5 Myr uncertainty in age, which propagates into a mass uncertainty.

Intensity ratios are the most rudimentary diagnostic of CO chemical depletion. Williams & Best (2014), Miotello et al. (2014), and Miotello et al. (2016) show how isotope-selective photodissociation, varying disk radii, and CO freezeout (or lack thereof) can affect intensity ratios, leading to a wide dispersion of C18O/13CO for disks of the same mass.

6.2 How line profiles reveal CO chemical depletion

The line profiles contain information not available in the integrated intensities. Here we demonstrate how a comparison of normalized line profiles from 13CO and C18O or C17O can diagnose CO chemical depletion, even over a factor of two range in disk mass. Figure 4 shows that the normalized line profiles of the rarer isotopologues become broader than those of the more common isotopologues as time proceeds. Since CO chemical depletion happens primarily in the outer part of the disk where Keplerian velocities are small, the fraction of radiation in the high-velocity line wings increases with time for optically thin lines. Optically thick lines can mask CO chemical depletion: the CO column density has to first decline to the $\tau\approx 1$ threshold (where $\tau$ is the optical depth in the line center) before intensity changes start to track abundance changes. For the $0.015\,{M_{\odot}}$ disk, the optically thin C17O and C18O emission tracks the changes in CO column densities well, so that the line profile gets significantly broader after $1$ Myr. 13CO ($J=3\rightarrow 2$) emission is optically thick well beyond the CO chemical depletion radius of $\sim 20\,{\rm AU}$, so does not reveal reductions in CO column density with time. Similar trends are found for the $0.03\,{M_{\odot}}$ disk. Because of the relatively high column densities, the change in C18O line profile is not obvious until $3\,{\rm Myr}$, but C17O still tracks the reductions in CO column density very well. Clearly, line profiles of rare isotopologues can diagnose CO chemical depletion .

However, emission line profiles reflect not only the CO/H2 abundance ratio, but the temperature and density structure of the disk as well. The temperature of our model disk decreases with time due to the young star’s dimming as it moves down the Hayashi track, and the density structure changes as the disk viscously evolves. To isolate the effect of CO depletion from density and temperature effects, we set up control models with a constant CO/H2 ratio throughout the disk and compare the emission line profiles produced by these with those produced by our fiducial, CO-depleted models. In the constant CO/H2 models, all atomic carbon available for gas-phase reactions is assumed to be in CO, and relative abundances of CO isotopologues are determined only by the abundance ratios of the isotopes (see Table 2 of Paper 1). The abundances normalized to the total proton density are $7.21\times 10^{-5}$ for CO, $9.34\times 10^{-7}$ for 13CO, $1.44\times 10^{-7}$ for C18O, and $3.13\times 10^{-8}$ for C17O.

Figure 9 shows side-by-side normalized line profiles for the $0.015$ M⊙ disk with evolving chemistry and the one with constant $f_{\rm CO}$. By 2 Myr, the profiles of rare isotopologues are noticeably wider than those of more common isotopologues and by 3 Myr, they are quite distinctive. These line profiles comparisons would provide clear evidence for ongoing CO chemical depletion and a potential way to correct for it (though there may be multiple rings of gas in some disks, which would complicate the line profile analysis; see Cleeves 2016). However, the actual (un-normalized) lines are very weak (Fig. 4), so diagnosing CO chemical depletion by comparing line profiles would be an expensive method in terms of observing time. Furthermore, variation in turbulent speeds between different layers of the disk could also produce different line profiles for 13CO, C18O, and C17O (Flaherty et al., 2015; Simon et al., 2015), an effect we have not explored here.

6.3 Spatial distributions of CO isotopologues

Our chemical evolution models also predict a characteristic radial dependence of the CO abundance. To translate this dependence into observables, we average the velocity-integrated, continuum-subtracted emission from our model disks in rings of $2$ AU in radius, and compare the spatial distributions of the $J=2\rightarrow 1$ emission line from the fiducial models to that from the constant CO models in Figure 10. The models shown in Figure 10 were run in face-on geometry for simplicity; all other simulated emission in the paper is calculated for $30^{\circ}$ inclination. The optical depth effects are enhanced in the face-on disks because the lines are not spread out by rotation. If a disk inclination angle is known, models could be run for that situation.

Results from the fiducial $0.015\,{M_{\odot}}$ disk are presented in the upper-left (early time), upper-right (2 Myr), and lower-left (3 Myr) panels. At the beginning of the evolution, the integrated intensities are lower in the fiducial model but the profiles are very similar to those in the constant CO disk for all isotopes. This is because part of the carbon is locked in CO2 in the fiducial model, whereas the constant CO model has all available carbon in CO. As the disk evolves, the intensities beyond 20 AU decrease dramatically in the fiducial model due to the depletion of CO, and the differences between the fiducial model and the constant CO model are the greatest for C17O lines, which have the lowest optical depth. We see sharp drops of intensity in models with $0.03\,{M_{\odot}}$ (lower right) as well. However because the chemical depletion happens more slowly in the $0.03\,{M_{\odot}}$ disk, the intensity profiles at 3 Myr of the $0.03\,{M_{\odot}}$ disk resemble those of the $0.015\,{M_{\odot}}$ disk at 2 Myr.

The effect of age is detectable with ALMA given enough integration time. We plot the time evolution of 13CO and C18O lines for the 0.015\mbox{M{}_{\odot}} disk in Figure 11, and the detection limit as a horizontal line at log $=-6.27$. At resolution, we get to 540 mJy in 1 km$/$s resolution in 10 hours. At $0\farcs 03$ spatial resolution, with a velocity channel of width 1 km$/$s , we reach a sensitivity of 540 mJy per beam in 10 hours. The line width for a perfectly face-on disk could be small, but it only takes a small inclination angle to significantly broaden the lines (for example, it takes only 6 degrees for Keplerian velocity to broaden the line to 1 km/s at 10 AU around a solar mass star). $0\farcs 03$ will resolve 4.2 AU at 140 AU and matches our resolution in the figure. For the purpose of illustration, we use the $J=2\rightarrow 1$ lines from a disk of $30^{\circ}$ inclination. The apparent edge of the disk in 13CO moves from beyond 70 AU to about 55 AU, and from beyond 70 AU to about 20 AU in C18O emission, well within the condensation front of CO in both cases. The migration of the “fake” snowline is a combined effect of the drop of CO emission intensities from the outer disk and the detection limit of the observation.

In disks where spatially resolved imaging is possible, a comparison between the C18O and the CO or 13CO spatial distribution should reveal CO chemical depletion, and may even provide enough information to re-construct the CO distribution with the help of a chemical model. Finally, after 3 Myr of evolution, we see that two rings of C2H have formed in the surface layers of the disk: a narrow ring at $1-3$ AU and a broader ring at $10-20$ AU. Bergin et al. (2016) found that C2H formation is possible only with a strong UV field and C/O$>1$, conditions which are replicated in our model disk surface layers due in part to the breakdown of CO molecules. We will explore the detectability of the C2H rings and the degree to which they may indicate CO depletion in future work.

7 Consequences for Observations

Our main result is that current observations and interpretations of CO isotopologues toward disks around solar-mass stars likely underestimate gas masses by substantial amounts. The surveys so far have emphasized shallow observations of many disks with modest spatial resolution. For example, Ansdell et al. (2016) surveyed $89$ protoplanetary disks in Lupus with ALMA and found very low gas masses when assuming an ISM CO abundance. The analysis of those observations has led to the conclusion that almost no disks have masses that exceed the MMSN. This result would appear to conflict with the high fraction of stars with evidence for planetary systems.

Our predicted line intensities are weak, and deriving correct masses requires deeper integrations to obtain either line profiles of weak lines or very high spatial resolution. Until those observations are available, we can offer only rough estimates of how much more massive the disks may be. As an example, we ask the following question: if we apply a simple correction for $f_{\rm CO}$, how many of the disks in Table 3 of Ansdell et al. (2016) might contain the mass of $0.01$ M⊙, the MMSN. If the mean age of stars in Lupus is $3\pm 2$ Myr, we should use \mbox{f_{\rm CO}}=0.136, the value at 3 Myr for the 0.015 M⊙disk, the closest to a MMSN. After applying this correction factor, 7 of 36 best guess masses exceed the MMSN, and 26 of 36 maximum masses exceed the MMSN, where best guess and maximum masses are defined by Ansdell et al. Within the uncertainties, a significant fraction of the disks in Lupus could contain enough mass to make a planetary system like ours. Recent analysis of similar observations toward disks in Cha I indicate a deficit of gas mass similar to, or worse than, that seen in the Lupus disks (Feng et al. 2017). For an age of 2 Myr and the 0.015 M⊙ disk, the value of $f_{\rm CO}$ is 0.18, which could increase the masses by a factor of about 5, leaving almost all disks in Cha I still short of the MMSN. Almost all these stars are lower in mass than our model star; models tuned to these stellar parameters would be needed to draw further conclusions.

CO chemical depletion also affects estimates of gas/dust mass ratios. Assuming an ISM-like CO/H2 abundance ratio, Ansdell et al. (2016) found a wide range of gas/dust ratios, with a median ratio of $\sim 15$, which is much smaller than the canonical value of $100$ observed in the ISM. Miotello et al. (2017) have shown that this problem persists when their models with radiative transfer and isotope-selective photodissociation are used.

We argue that the low gas-to-dust ratios measured in Ansdell et al. (2016) are likely a result of a low CO/H2 ratio due to CO chemical depletion. In our 0.015 M⊙disk model, the closest to a MMSN, only 13.6$\%$ of carbon is contained in CO gas by 3 Myr. If the mean age of stars in Lupus is $3\pm 2$ Myr, and we apply the simple correction for CO fraction ($f_{\rm CO}$= 0.136) to the Lupus data with both 13CO and C18O detections, we obtain gas-to-dust ratios of 53 to 1700, with a median gas-to-dust ratio of 108. The gas-to-dust ratios correcting for the fraction of C in CO in our models are plotted in Fig.  12.

The high values of the gas-to-dust ratios in the $f_{\rm CO}$ corrected gas-to-dust ratios could be a combined result from grain growth, weak-ionizing environments, and young disk ages. Gas/dust mass ratios are not necessarily the same as gas/solid mass ratios: the apparent dust mass is typically calculated from sub-mm dust emission, which is mostly sensitive to dust grains around sub-mm size or smaller. Once dust grains grow into cm-size pebbles, the reduction of observable dust mass could drive the gas/dust mass ratio to larger than the canonical molecular cloud value of $100$. If only $10\%$ of the solid mass is in the form of dust observable in $890$ µm, as we assume in our models, one would have to correct for the “dust fraction” ($f_{\rm dust}$) when measuring the gas-to-dust ratios. We apply the correction for dust fraction on the $f_{\rm CO}$ corrected data, and plot the gas-to-dust ratios as blue squares in Fig.  12. Most the disks with high gas-to-dust ratios from the $f_{\rm CO}$ corrections can be come consistent with the ISM value of 100 if additional corrections for grain growth are included.

Instead of assuming that dust evolution or disk clearing drives the variation in observed Lupus gas/dust ratios, we point to the large age spread of pre-main-sequence stars in Lupus. Ansdell et al. quote 1-3 Myr, while Mortier et al. (2011) find either 0.1 to $>15$ Myr or 0.3 to $>15$ Myr depending on the model isochrones (though the $>10$ Myr objects were not detected by Ansdell et al. and are therefore not represented in Figure 12). The large spread of gas/dust ratios calculated by Ansdell et al. is likely the result of the variation of CO/H2 abundance ratio with protostellar age. Given the uncertainties in both gas and solid mass, the data so far do not rule out ISM gas/dust ratios.

There are a few caveats to bear in mind regarding these simple mass correction factors. First, real disks evolve at different rates: a disk with a high incident cosmic-ray flux loses CO more quickly than a disk shielded by a “T-Tauriosphere” that deflects cosmic rays (Cleeves et al., 2013, 2015). Even if star ages could be measured perfectly, the ages alone do not give enough information about the chemical evolutionary stage of the disk to compute a CO/H2 ratio. In other words, our line profile diagnostics are measuring a “chemical age.” Comparison of 13CO and C18O line profiles gives an empirical diagnostic of the level of CO chemical depletion .

Second, chemical depletion of CO operates on a million-year timescale. If episodic accretion of the type observed in FU Orionis outbursts happens during the T-Tauri phase at intervals smaller than $\sim 1\,{\rm Myr}$ (e.g., Dunham et al., 2010; Kim et al., 2012; Martin et al., 2012; Green et al., 2013), desorption of CO2 ice followed by CO2 dissociation from cosmic ray-induced photons could raise the CO gas abundance throughout the disk, not just in the inner 20 AU as in our current model disk. Cieza et al. (2016) have already proposed chemical alteration in the disk surrounding V883 Ori, which is now mid-outburst. We have not modeled FU Orionis outbursts—our model T-Tauri star evolves smoothly on the Hayashi track.

Other processes may further deplete CO in the outer regions of disks. For example, Xu et al. (2017) describe a process of “runaway freeze-out”, which considers vertical transport at fixed radii. They find that higher layers, too warm for freeze-out, can become depleted by transport to the colder, lower layers. Kama et al. (2016a) suggest the same mechanism to explain the low atomic C and O abundances in the TW Hya disk. Applied to CO, the vertical transport/freezeout effect would primarily act at larger radii than we consider and affect primarily the emission from the more common, hence more optically thick, isotopes. Applied to the complex organics that sequester the carbon in our models, runaway freezeout could remove even more carbon from the gas phase, allowing for even higher H2 masses to be consistent with observed gas-phase carbon abundances.

Our model’s growing CO abundance in the inner disk and CO chemical depletion in the outer disk would be observationally almost identical to a disk with CO frozen onto grain surfaces in the outer disk, followed by inward radial drift of the grains and desorption of CO in the warm inner disk. The radial-transport effect, first applied to water ice, has been suggested on theoretical (Ciesla & Cuzzi, 2006; Du et al., 2015) and observational (Hogerheijde et al., 2011; Kama et al., 2016b) grounds, diagnosed from the presence of rings of small hydrocarbon chains (Bergin et al., 2016) and ammonia gas (Salinas et al., 2016) in the outer disk of TW Hya, and inferred in cases where C/O$>1$ (Bergin et al., 2016). Robust temperature measurements in the outer disk would be required to distinguish between chemical CO chemical depletion and CO freezeout followed by grain inspiral.

8 Conclusions

Our key findings and suggestions for observing strategies are summarized below.

CO abundance varies both with distance from the star and as a function of time as CO is dissociated and the carbon gets sequestered in organic molecules that freeze onto grain surfaces (chemical depletion) on a million year time scale (§2). CO chemical depletion will cause very large underestimates in gas mass and gas-to-dust ratios when CO observations are used. One would need to correct for the chemical depletion of CO in order to correctly estimate the disk gas mass. The CO abundance correction factor ranges from 3 to 8 for the models we have run (§4). 2. 2.

Even though CO is destroyed by ionized helium throughout most of the disk, it has a higher-than-interstellar abundance in the inner 20 AU of the disk, where the temperature is relatively high. Adopting a constant $T=20$ K for the CO reservoir will underestimate gas masses by a further factor of about 2 (§4). 3. 3.

The high CO abundance in the inner disk also results in high optical depth even for the mostly optically thin isotopologue we investigated (C17O), and one could underestimate the disk gas mass due to the optical depth effect even after correcting for the CO chemical depletion in the outer disk and using the correct CO-averaged temperature (§4). 4. 4.

The CO-abundance time evolution also introduces a disk age-mass degeneracy—a massive, old disk may look similar to a young, less massive one in CO rare isotopologue emission (§5). 5. 5.

One can diagnose CO chemical depletion by comparing the line intensity ratios (§6.1) or emission line profiles of multiple isotopologues (§6.2). If the disk is spatially resolved, one can also use spatial distribution of CO beyond $20$ AU to probe the chemical depletion of CO. (§6.3). 6. 6.

Very different CO optical depths in different parts of the disk produce a complicated rotation diagram that does not probe the disk temperature well. One would underestimate the average temperature of CO molecules by deriving it from the lowest two J values (§C.1). 7. 7.

Higher-$J$ lines underestimate the disk mass by more than do lower-$J$ lines. The higher-$J$ lines miss low temperature gas. We suggest using low excitation lines (e.g., $J=1\rightarrow 0$) to estimate the disk mass to minimize the temperature and optical depth effects. (§C.2). 8. 8.

The strategy that comes closest to recovering the correct mass for our models is to use the formulae from Miotello et al. (2016) and to divide by the $f_{\rm CO}$ from our models (§4). 9. 9.

If we correct the “best-guess” gas masses in Ansdell et al. (2016) by the smallest value of $f_{\rm CO}$ (the value at 3 Myr for the 0.015 M⊙ disk), 7 of the disks could have masses of the MMSN; if we apply our $f_{\rm CO}$ correction to the maximum masses, 26 of 36 could reach the MMSN (§7). The age of Lupus is $3\pm 2$ Myr (Alcalá et al., 2014), so this correction is reasonable. A recent survey of disks in Cha I suggest almost no disks with MMSN masses. While still problematic for planet formation models, correcting for CO chemical depletion suggests better numbers than finding that no disks contain the MMSN. 10. 10.

Given reasonable uncertainties in both gas and solid masses, many observed disks can have gas/dust ratios consistent with the ISM value of 100 (§7).

Work by MY, KW, SDR and NJT was supported by NASA grant NNX10AH28G and further work by MY and SDR was supported by NSF grant 1055910. This work was performed in part at the Jet Propulsion Laboratory, California Institute of Technology. NJT was supported by grant 13-OSS13-0114 from the NASA Origins of Solar Systems program. MY was supported by a Continuing Fellowship from the University of Texas at Austin. We acknowledge helpful input from L. Cleeves, E. Bergin, E. F. van Dishoeck, K. Öberg, J. Huang, M. Ansdell, and A. Kraus. We are grateful for Feng Long et al. sharing their paper in advance of publication, and Huang et al. sharing processed ALMA data for our model comparison.

Appendix A Sensitivity to model parameters

The effect of our assumption of LTE in the CO level populations and the fixed vertical grid are explored in the next two sub-sections.

A.1 LTE vs. NLTE

LIME is capable of computing energy level populations either in local thermodynamic equilibrium (LTE) or in the more complex non-LTE case, where the kinetic temperature and excitation temperature are different. In Paper 1, the lowest density found in the modeled region within 70 AU of the star is about $10^{9}$ hydrogen molecules per cubic centimeter. Even in the most diffuse regions of our model disk, the density should be high enough for collisions to dominate CO excitation so that LTE is a good approximation for the energy level population.

We demonstrate that the LTE approximation is appropriate for our model disk by comparing the results of LTE and non-LTE models for different CO isotopologues at multiple epochs. For each emission line of each isotopologue, the line profiles from the two models are indistinguishable. In figure 13 we show the worst case example of C17O J=$6-5$ emission at the beginning of the disk evolution. The high excitation energy of the J=$6-5$ transition means the LTE approximation requires a high collision rate and is hardest to satisfy, yet the differences between the LTE and non-LTE line profiles are still negligible.

A.2 Definition of the disk surface

Another parameter that can influence the computed line profiles is the placement of the disk surface. In the disk models of Landry et al. (2013), on which the chemical models in Paper 1 are based, the top surface of the computational grid is placed at the layer where the Rosseland mean optical depth to the disk’s own radiation, integrated downward from infinity, is 0.2. The $\tau=0.2$ surface is located between one and two pressure scale heights above the disk midplane, depending on the distance from the star and time. Since most of the low-J CO emission comes from the disk interior where the bulk of the disk mass is located, the height of the grid surface should have minimal effect on the line profiles. For higher values of J, more of the emission would come from the warm surface and modeling line profiles accurately would require the computational grid to extend well into the tenuous disk atomsphere. Here we assess the effects of grid surface placement on our model emission line intensities and profiles.

The chemical model in Paper 1 uses a static grid that does not evolve with time. Grid cells are spaced logarithmically in radius $r$ and linearly in aspect ratio $z/r$. At the beginning of evolution, which roughly corresponds to the start of the T-Tauri phase, the disk has a high scale height due to heating from the luminous protostar. The scale height decreases throughout the $3\,{\rm Myr}$ evolution due to both protostellar dimming and viscous dissipation. As the disk flattens, grid layers with high $z/r$ begin to empty out. To test the effect of surface placement on the radiative transfer calculation, we construct a new disk model by artificially removing the top three filled (not empty) $z/R$ layers from each time snapshot of our fiducial model from Paper 1. At six different epochs, we compute line profiles for the $1\rightarrow 0$, $2\rightarrow 1$, $3\rightarrow 2$, $4\rightarrow 3$, $5\rightarrow 4$, and $6\rightarrow 5$ transitions for all CO isotopologues using both the new “remove-top” model and the fiducial model. Removing emitting layers affects the peak intensities of the optically thin C17O emission lines more than any other isotopologue. To verify that our radiative transfer models capture essentially all of the disk emission, we demonstrate how the grid surface placement affects C17O.

At the beginning of disk evolution, the top three grid layers contain very little mass and have little effect on the emission line intensity or profile (Figure 14, left panel). As the the disk cools and the top grid layers empty out, the mass contained in the filled grid layers increases. Line profiles from the fiducial and remove-top models differ the most after 3 Myr of evolution (Figure 14, right panel). Yet even at 3 Myr, the difference in the $J=3\rightarrow 2$ peak intensity between the remove-top and the fiducial model would be difficult to distinguish in observations, and the normalized line profiles are almost identical.

The differences in peak intensities between the fiducial and remove-top models become significant for higher-J emission due to the higher energy needed to populate the upper state. For example, the peak intensity for $J=6\rightarrow 5$ differs by roughly $33\%$ between the fiducial and remove-top models at 3 Myr of evolution. However, lower J lines ($J=1\rightarrow 0$ to $3\rightarrow 2$) are more commonly used for disk studies. In this paper, we focus our line profile discussion on $J=3\rightarrow 2$ and $J=2\rightarrow 1$ transitions, which are the most observationally relevant transitions. For our purposes, we have adequately modeled the CO rare isotopologue emission, even though our models do not extend vertically to a large number of scale heights.

Appendix B Measuring disk mass from the integrated intensities of CO isotopologues

For optically thin emission lines, the total number of molecules per degenerate sublevel in the upper state can be written simply as

[TABLE]

where $L_{J}$ is the line luminosity, $g_{J}$ is the degeneracy, $A_{J}$ is the Einstein A value for the transition from upper state with quantum number $J$, $h$ is Planck’s constant, and $\nu$ is the frequency of the transition. $L_{J}$ is related to the integrated line flux $F_{l}$ as

[TABLE]

where $D_{\rm cm}$ and $D_{\rm pc}$ are the star’s distance in cm and pc, respectively, and $c$ is the speed of light. For the last part of Eq. B2 and Eq. B4, Fl is in the unit of ${\rm Jy-km\ s^{-1}}$, Dpc is in the unit of pc, and all the other values are in cgs units. The Einstein A value for rotational transitions from level $J$ to level $J-1$ is given by

[TABLE]

where $|\mu(J,J-1)|$ is the electric dipole matrix element and $\mu_{\rm D}$ is the dipole moment measured in Debye ($10^{-18}$ esu-cm). Combining these, we can write

[TABLE]

where $B_{\rm GHz}$ is the rotation constant in GHz. The values of $B_{\rm GHz}$ vary slightly with the isotopologue and are easily obtained from on-line sources, but for reference, \mbox{B_{\rm GHz}}\approx 55 to 58 GHz for the isotopes discussed here.

Appendix C Does a single temperature characterize the CO emission?

C.1 Temperature Estimates from Rotation Diagrams

A temperature estimate is required in order to evaluate the partition function that will translate observed intensities into gas column densities. While the disk temperature can reach $\sim 1400$ K at the dust sublimation front (Muzerolle et al., 2003), Andrews & Williams (2005) and Andrews et al. (2013) suggest that most submillimeter emission comes from dust at $\sim 20-25$ K and a fixed temperature of 20 K for the gas is often used in the simplest methods to measure mass. Here we will use plots of level populations versus energy of the level above ground—called rotation diagrams for rotational transitions—to test whether CO in our model disk is well characterized by a single temperature. We calculate the rotation diagrams for the $0.015\,{M_{\odot}}$ disk and the $0.03\,{M_{\odot}}$ disk using the equations in Appendix B, following the method of Green et al. (2013). Note that the emission must be optically thin (which is often assumed for CO isotopologues) for luminosity to be directly proportional to the upper state population.

Figure 15 shows the 13CO rotation diagram of our $0.015\,{M_{\odot}}$ model disk (left, circles) and our $0.03\,{M_{\odot}}$ model disk (right), both after 2 Myr of evolution. A gas reservoir with a single temperature would yield a straight line in $\ln(N_{J}/g_{J})$ as a function of $E_{\rm up}$. Unfortunately, Figure 15 shows that no single temperature characterizes the CO level populations, consistent with Herschel observations of Herbig Ae/Be and T-Tauri disks by Meeus et al. (2013), Fedele et al. (2013), van der Wiel et al. (2014) and Fedele et al. (2016). Even when we construct an artificial disk with the same density and abundance structure as our 2 Myr models, but with a constant temperature enforced at all points (Figure 15, squares and triangles), the rotation diagram still appears to come from a disk with a range of temperatures because of the varying optical depth of different transitions. If we derive a temperature from the lowest two J values, the temperature is quite low—about 10 K. CO rotation diagrams, even for rare isotopologues, are not reliable ways to measure protostellar disk temperatures.

C.2 Inferred Mass Depends on the Transition Used

Because the temperature and optical depth are different in various parts of the disk, the mass estimation also depends on which isotopologue and which emission line is used. We show the mass estimated by different transitions and isotopologues at $100$ yr, $2$ Myr, and $3$ Myr of the disk evolution in Figure 16. All models assume optically thin emission for simplicity. The estimated mass decreases as we use higher-J lines with higher excitation energy for the mass estimation. This is partially contributed by the large optical depth in the inner hot regions of the disk, partially because we are missing the low temperature CO that does not emit much at higher-J. Our models suggest that it is best to use lower $J$ transitions to measure disk mass.