Mass inventory of the giant-planet formation zone in a solar nebula analog

TL;DR

This study provides the first spatially resolved measurements of bulk gas mass distribution in a protoplanetary disk, offering new insights into planet formation zones and dust growth processes.

Contribution

It introduces a novel method combining $^{13}$C$^{18}$O and C$^{18}$O observations to directly measure gas mass distribution in a protoplanetary disk.

Findings

Gas mass distribution in the 5-21 AU zone of TW Hya was quantified.

The gas-to-dust mass ratio was found to be approximately 140.

The radial gas distribution aligns with a viscous disk model but is flatter than in our Solar System.

Abstract

The initial mass distribution in the solar nebula is a critical input to planet formation models that seek to reproduce today's Solar System. Traditionally, constraints on the gas mass distribution are derived from observations of the dust emission from disks, but this approach suffers from large uncertainties in grain growth and gas-to-dust ratio. On the other hand, previous observations of gas tracers only probe surface layers above the bulk mass reservoir. Here we present the first partially spatially resolved observations of the $^{13}$C$^{18}$O J=3-2 line emission in the closest protoplanetary disk, TW Hya, a gas tracer that probes the bulk mass distribution. Combining it with the C$^{18}$O J=3-2 emission and the previously detected HD J=1-0 flux, we directly constrain the mid-plane temperature and optical depths of gas and dust emission. We report a gas mass distribution of 13$^{+8}_{-5}\times$(R/20.5AU)$^{-0.9^{+0.4}_{-0.3}}$ g cm$^{-2}$ in the expected formation zone of gas and ice giants (5-21AU). We find the total gas/millimeter-sized dust mass ratio is 140 in this region, suggesting that at least 2.4M_earth of dust aggregates have grown to >centimeter sizes (and perhaps much larger). The radial distribution of gas mass is consistent with a self-similar viscous disk profile but much flatter than the posterior extrapolation of mass distribution in our own and extrasolar planetary systems.

Tables

Table 1. Table 1: Best fitting TW Hya disk model parameters

Rsnow	Tsnow	$q_t$	${\tau }_0$	$q_1$	$q_2$	log$\delta$	12C/13C	Rring	Aring	${\sigma }_{\mathrm{ring}}$
(AU)	(K)		(km s-1)					(AU)		(AU)
20.5$\pm 1.3$	27${}_{-2}{}^{}{}_{}{}^{+3}$	-0.47${}_{-0.07}{}^{}{}_{}{}^{+0.06}$	2.5${}_{-0.5}{}^{}{}_{}{}^{+0.6}$	-0.18${}_{-0.23}{}^{}{}_{}{}^{+0.39}$	-1.9$\pm 0.1$	-0.94$\pm 0.14$	40${}_{-6}{}^{}{}_{}{}^{+9}$	60${}_{-1}{}^{}{}_{}{}^{+1}$	6.3${}_{-0.8}{}^{}{}_{}{}^{+1.2}$	11$\pm 1$

Table 2. Supplementary Table 2: Fixed model parameters for the disk of TW Hya

Stellar mass	Distance	Inner radius	Outer radius	inclination		Turbulence	Position angle
M${}_{}{}^8{}_{\star }{}^{}$	$d$[7]	rin[48]	rout	$i$[49]	VLSR[50]	$v_{\mathrm{turb}}$	PA[50]
(M⊙)	(pc)	(AU)	(AU)	(deg)	(km s-1)	(km s-1)	(deg)
0.6	55	2	150	7	2.8	0.01	155

Equations

$$\{\begin{array}{rl}{I}_d& =B_{\nu }(T)(1-e^{-{\tau }_d})\\ I_{{\mathrm{C}}^{18}\mathrm{O}}& =B_{\nu }(T)(1-e^{-({\tau }_d+{\tau }_{C^{18}O})})-I_d\\ I_{{}^{13}{\mathrm{C}}^{18}\mathrm{O}}& =B_{\nu }(T)(1-e^{-({\tau }_d+{\tau }_{{}^{13}{C}^{18}O})})-I_d\\ & =B_{\nu }(T)(1-e^{-({\tau }_d+{\tau }_{C^{18}O}/ratio)})-I_d\end{array}$$

Full text

{affiliations}

Department of Astronomy, University of Michigan, 1085 S University Ave., Ann Arbor, MI 48109, USA

Division of Geological & Planetary Sciences, MC 150-21, California Institute of Technology, 1200 E California Blvd., Pasadena, CA 91125, USA

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

The primary theory for the formation of giant planets is the so-called core accretion scenario, where a rock+ice core forms through the coagulation of planetesimals until it becomes sufficiently massive to accrete a gaseous envelope[1]. In this theory, the spatial distribution of gas in the primitive nebula is not only critical to the later accretion of the atmosphere of giant planets, but also plays an important role in the early planetesimal formation processes as the transport and mixing of grain aggregates depend on the gas turbulence and density[5, 6].

Here we report a high spatial and spectral resolution observational study of the gas mass distribution in the expected formation zone of gas and ice giants using the Atacama Large Millimeter/submillimeter Array (ALMA). In particular, we robustly detect the 13C18O and C18O $J$=3-2 line emission (5$\sigma$ and 29$\sigma$, respectively) in the closest protoplanetary disk around the T Tauri star TW Hya (see Figure 1). TW Hya, at a distance of only$\sim$55 pc[7], is one of the best studied disk systems and makes a useful analog to the solar nebula: the central star has a stellar mass about 0.6-0.8 M⊙[8]; the disk itself has a total mass (gas and dust) of $\sim$0.05 M⊙[9], comparable to the total mass of 0.01 M⊙ in the Minimum Mass Solar Nebula (MMSN) which is the lowest possible mass needed to reproduce the Solar System. The age of TW Hya is 8$\pm$4 Myr[10], within the expected window for giant planet formation[11].

These observations were specifically designed to probe the mid-plane interior to 30 AU. CO, one of the most volatile species after H2 and He, remains in gas phase even in the mid-plane within a large radius[12]. Previous observational studies have used more abundant CO isotopologues (e.g. 13CO), which have stronger line emission, to study the gas mass distribution in protoplanetary disks[4]. However, inside the mid-plane CO snowline, more abundant CO isotopologue line emission primarily arises from layers at much larger scale heights. As a result, the derived gas mass distribution is highly model dependent, hinging on the adopted temperature and vertical density structures. In Figure 2, we show the vertical emitting regions of CO $J$=3-2 line emission from four isotopologues under representative properties for the TW Hya disk. It shows that only the 13C18O (3-2) line emission arises from both sides of the disk and thus provides a robust probe of the bulk mass.

To accurately constrain physical properties in the mid-plane, we employ a parameterized disk model to reproduce the observations. A full description of the model and the Markov chain Monte Carlo (MCMC) fitting algorithm is given in the Methods section. The critical assumption of the model is that the 13C18O and C18O  emission arises predominately from the mid-plane or immediately adjacent regions, and therefore that the emitting gas can be characterized by a one-dimensional mid-plane temperature. This assumption is reasonable since C18O and 13C18O are 10*-10*–10*-6* times less abundant than H2[13], and the disk is nearly isothermal within one scale height[14]. To include potential extinction from dust, we treat the C18O, 13C18O line emission and dust continuum simultaneously. For the dust emission, we employ the radial profile of surface brightness derived by Andrews et al. (2016)[15], because these continuum observations were taken at a similar frequency as our continuum data (345 v.s. 320 GHz) but with much higher spatial resolution ($\sim$1 AU). We use simple power-law functions to characterize the radial profiles of the temperature and integrated line strength. The best-fitting parameters and their uncertainty are listed in Table 1, and the best-fitting model is plotted in Figure 1 in comparison with the data.

We report the first direct measurement of the mid-plane CO snowline at 20.5$\pm 1.3$ AU, in the TW Hya disk. The 13C18O $J$=3-2 emission is optically thin ($\tau<$0.2) throughout the disk, and the dust continuum emission has $\tau\sim$0.5, suggesting that the 13C18O line emission arises from both sides of the disk and thus provides model-independent constraints on the mass distribution. The C18O emission is optically thick ($\tau\sim$8) inside 20.5 AU and therefore provides direct constraints on the gas temperature rather than the mass distribution (Supplementary Figure 4). The observations trace the mid-plane between 5-20.5 AU: the outer limit is the mid-plane CO snowline location (beyond that CO condenses at the mid-plane); the inner boundary is limited by the radius of 10% cumulative flux at the highest velocity channel with a 13C18O  emission detection (Methods and Supplementary Figure 5).

With the temperature and optical depths from the best-fitting model, we can derive the column density distribution of 13C18O and then the gas mass distribution for a given 13C18O abundance ($x_{\rm{}^{13}C^{18}O}$). Here we employ the spatially unresolved HD $J$=1-0 line flux of the TW Hya disk to estimate the $x_{\rm{}^{13}C^{18}O}$[9]. Combining the uncertainties in our model parameters and in the HD line flux, we find $x_{\rm{}^{13}C^{18}O}$=1.7${}^{+1.3}_{-0.8}\times$10*-10* and $\Sigma_{\rm gas}$=13${}^{+8}_{-5}\times$(R/20.5 AU)${}^{-0.87^{+0.38}_{-0.26}}$ g cm*-2* between 5-20.5 AU (for details, see Methods).

For dust mass, the conversion from optical depth to surface density depends on the dust opacity $\kappa_{\nu}$, which may vary with radius due to different levels of dust growth. The growth level of dust aggregation is usually estimated with a parameter $\beta$ where $\kappa_{\nu}\propto\nu^{\beta}$. Combining our best-fitting temperature structure and previous continuum observations at 145, 233, and 345 GHz[15, 20], we find that the $\beta_{\rm 145-233\,GHz}$ increases from $\sim$ 0 to 1.1 between 5-20.5 AU, consistent with the finding of Tsukagoshi et al. (2016) which used a different temperature structure. But we also find that the $\beta_{\rm 233-345\,GHz}$ is roughly flat in the same region. As the dust emission in this region has $\tau<1$ (Supplementary Figure 7), the $\beta\lesssim 1$ suggests that significant fraction of dust has grown to mm-size or larger[21], and the discrepancy of $\beta$ among different wavelengths hints that the dust size distribution in this region is possibly significantly different from the typically assumed power-law function with an index of -3.5. Because detailed modeling of the dust size distribution is beyond the scope of this paper, here we adopt a constant $\kappa$(0.93 mm) = 3.1 cm2g*-1* to estimate the dust mass of grains no larger than millimeter. The opacities adopted are similar to those used in previous models. Assuming a constant $\kappa_{\nu}$, the radial distribution of mm-sized dust mass between 5-20.5 AU is shown in Figure 3.

We measure the total gas-to-mm-sized dust ratio is 140${}^{+43}_{-31}$ for the 5-20.5 AU region, much higher than the typical ratio of 100 in molecular clouds. To recover a gas-to-dust ratio of 100, we find that 2.4${}^{+3.2}_{-2.0}$ M*⊕* solid masses are missing. This estimation is actually a lower limit, because the dust radial drift to the inner region should in general decrease the total gas-to-dust ratio in the inner disk[6] and because the overall gas surface density also decreases with time. Therefore more than 2.4 M*⊕* masses of solids may have grown to at least centimeter-sized and perhaps much larger bodies in the TW Hya disk. This estimation sheds new light on the previous skepticism of planet formation in the TW Hya disk. Recently, a series of dark rings were discovered in the (sub)mm continuum emission of the TW Hya[15] and other disks[22], which has excited heated debates on the origin of these mysterious rings and especially whether they formed via planet-disk interaction[23, 24]. Simulations suggest that the rings are too shallow and narrow to be gaps carved by Jupiter-like planets, but smaller bodies are still possible[25].

The direct measurement of gas mass distribution in the TW Hya disk provides the first data point to compare with the solar nebula over a similar radial region. Compared to the traditional MMSN, the TW Hya disk has less mass between 5-21 AU, though the central star is slightly less massive than the Sun. More importantly, the mass distribution in TW Hya is much flatter in this region (Figure 3). Extending the comparison to exoplanetary systems, the MMSN has its counterpart the MMEN, which is based on statistics of the close-in super-Earth population ($<$0.4 AU). The mass distribution in the MMEN models is at least as steep as the MMSN and perhaps much steeper (from -1.6 to -2.35), depending on the assumption of feeding zone size[18, 19].

Although more observations are needed to determine whether the gas mass distribution in TW Hya is common in protoplanetary disks, its shallow nature follows that derived from the dust emission[2, 3], and the apparent difference questions the use of MMSN/MMEN models as a standard initial mass distribution. The fundamental assumptions of these models are that planets formed at their current location and that the metal-to-(H+He) elemental ratio is solar throughout the solar nebula. In fact, the in situ formation assumption has been challenged by the discovery of hot Jupiters in extrasolar systems, where these close-orbit giant planets are too massive to form at their current location and probably are results of inward migration[26]. However, Solar system formation models that include planet migration require even steeper initial mass distribution than the MMSN[27]. If the gas mass distribution of the solar nebula was similar to the TW Hya disk, this large discrepancy means that the mass distribution of planetesimals (represented by planetary mass distribution before migration) must entirely decouple from the gas mass distribution in the solar nebula. In short, the MMSN/MMEN models may offer useful estimation as to the total mass of the natal protoplanetary disks but should not be used as quantitative templates of initial mass distribution for planet formation models.

Ultimately, the quantitive prediction of the mass distribution in protoplanetary disks should come from a detailed understanding of angular momentum transport in disks, for which observational studies such as this provide useful constraints. For example, in the self-similar solution of accretion disks[28], the power-law index $\gamma$ of the radial variation of viscosity (where $\nu\propto r^{\gamma}$) equals the power-law index ($p$) of surface density ($\Sigma_{\rm gas}\propto r^{-p}$ in the region r$\ll$ characteristic radius of the disk ). If the disk viscosity is mainly driven by turbulence (presumably due to magnetorotational instability), the effective kinematic viscosity is usually characterized by the Shakura-Sunyaev $\alpha$-model: $\nu=\alpha c_{s}^{2}/\Omega$, where $c_{s}$ is the isothermal sound speed ($c_{s}^{2}\propto T_{\rm gas}$) and $\Omega$ is the Keplerian angular velocity. Taking our measurement of $T_{\rm gas}$ and $\Sigma_{\rm gas}$, we find an $\alpha\propto r^{0.16\pm 0.39}$ between 5-21 AU. This result is consistent with the typically used constant $\alpha$ models, within uncertainties.

With the TW Hya measurement alone, however, it is not yet possible to test realistic magnetohydrodynamic models of disk evolution, because current simulations predict that the surface density profile varies widely depending on the input parameters and age of the disk[29, 30]. In the next few years, these models will be tested by measuring mass surface density profiles in a sample of disks with different ages and masses using techniques such as those describe here.

0.1 ALMA observations.

TW Hya was observed with the Atacama Large Millimeter/submillimeter Array (ALMA) as a part of the Cycle 3 proposal 2015.1.00308.S (P.I. Bergin). Observations were carried out on 2015 March 8, with 40 antennas and baselines of 15-460 m (16-492 k$\lambda$). The total on-source integration time was 67.5 min. The nearby quasar J1058+0133 was used for bandpass calibration, and J1037-2934 for both flux and phase calibration. The correlator was configured to simultaneously observe four spectral windows (SPWs). Three of the SPWs were used for line observations with a resolution of $\delta\nu$ = 122 kHz (a spectral resolution of 0.23 km s*-1* after the default Hanning smoothing is applied), and one for dedicated continuum observations. The 13C18O $J$=3-2 line was placed in SPW 1, and C18O 3-2 line in SPW 2 and 3. The three line SPWs are all 117.2 MHz wide and start at 314.171, 329.319 and 329.264 GHz, respectively. The SPW 0 was used for continuum observations with a total bandwidth of 1.875 GHz.

The visibility data were calibrated by ALMA staff following standard procedures with the pipeline. We further performed self-calibration on the calibrated data in phase and amplitude in CASA 4.6.12, using the dedicated continuum observations in SPW 0. The self-calibrated phase and gain solutions were then applied to the three line SPWs. The data were subsequently continuum subtracted using line free channels. We generated synthesized images using a CLEAN deconvolution with a manual mask (Briggs weighting with a robust parameter of 0.5). For the 321 GHz continuum emission map, the synthesized beam is 0.54${}^{\prime\prime}\times 0.47^{\prime\prime}$, with a position angle of -87 degrees. The C18O line image cube has a noise level of 6.5 mJy beam*-1* in each 0.23 km s*-1* wide velocity bin, and the 13C18O line image cube has a noise level of 4 mJy beam -1. The peak signal-to-noise ratio in the moment zero maps is 5 for 13C18O  and 29 for C18O, respectively. The synthesized images of line and continuum emission are displayed in Figure 1.

The observations of 13C18O and C18O were designed to probe the mid-plane interior to 30 AU. In the radial direction, CO remains in gas phase even in the mid-plane within a large radius[12]. In the vertical direction, previous studies have shown that the most abundant CO isotopologues (12CO and 13CO) readily become optically thick near the surface layer of protoplanetary disks and thus do not provide constraints on the deep interior [31, 32]. The two CO isotopologues chosen, 13C18O and C18O, are on average 38,433 and 577 times less abundant than 12CO in the interstellar medium (ISM)[13]. Their optical depths are correspondingly much smaller than the 12CO and 13CO lines, and thus can probe much deeper regions in protoplanetary disks.

0.2 Parameterized disk model.

Our goal is to constrain three key physical properties of the inner region of the TW Hya disk: the location of the mid-plane CO snowline (rsnow), the mid-plane temperature structure inside rsnow, and the radial mass distribution of the gas. To constrain these physical properties, we use an axisymmetric 1-Dimensional model to reproduce the 13C18O and C18O $J$=3-2 line emissions and dust continuum emissions at 345 GHz. The axisymmetric approximation is valid since no clearly asymmetric features have been seen in the line and continuum emissions from the TW Hya disk to date. The model structure consists of two distinct radial zones: inside and outside the mid-plane CO snowline (see Supplementary Figure 6). The two zones are used to model the break of the CO gas column density distribution due to condensation.

Our key focus is the nature of the region inside rsnow. The rsnow location is measured from the size of the bright central peaks in the 13C18O and C18O  images. We assume that inside rsnow, the C18O, 13C18O 3-2 line and continuum emissions primarily arise from the mid-plane or immediately adjoint regions, and therefore that the emitting gas column at each radius can be characterized by the mid-plane temperature. We will discuss the self-consistency of this assumption for our best-fitting model later. The abundance ratio of C18O/13C18O is assumed to be constant throughout the disk and modeled as a free parameter[33].

Under these two assumptions, the temperature and optical depths of line/continuum emission can be independently constrained by the 13C18O and C18O $J$=3-2 line and continuum emissions. The observed C18O line emission at the image center is only eight times stronger than the 13C18O line emission, much smaller than the typical ISM ratio of $\sim$69. The C18O line emission is thus likely to be optically thick, which provides direct constraints on the temperature inside the rsnow (see Supplementary Figure 4). With the temperature structure constrained, the optically thin 13C18O line emissions provide constraints on the CO gas column density distribution, and the continuum emissions constrains the dust optical depth. The CO gas column density distribution derived is then used as a proxy of the gas mass distribution, as current chemical models suggest that the CO/H2 abundance ratio is roughly constant in disk regions warmer than the CO condensation temperature[34]. The absolute scale factor is constrained by the previously observed HD (1-0) line flux[9].

The detailed construction of the model is described below. Inside rsnow, the continuum-subtracted line intensity and continuum intensity at each radius are modeled as an isothermal slab with the following equations:

[TABLE]

where $\tau_{\rm C^{18}O}(r,v)=\tau_{0}\times(r/r_{\rm snow})^{q_{1}}\times\phi_{v}$ (see Supplementary Figure 6), and $\phi_{v}$ is the local line profile that is modeled as a Gaussian with a standard deviation  $\sigma_{\rm V}=\sqrt{2k_{\rm B}T/m_{\rm CO}+v_{\rm turb}^{2}}$. We fix $v_{\rm turb}$ at 0.01 km s*-1*, a turbulence level consistent with an $\alpha$ of 10*-3* at 30 K, as previous ionization studies of the TW Hya disk[35] suggest that the turbulence level in the inner 50 AU region is likely to be very low. In addition, we test several models with different values of $v_{\rm turb}$. We find that using smaller values of $v_{\rm turb}$ does not change our best fitting results. Using a much larger $v_{\rm turb}$[36] (0.1 km s*-1*) would make the surface density profile slightly shallower but the difference is three times smaller than our uncertainty. We model the macroscopic velocity field as a Keplerian velocity field[37], and adopt stellar properties and disk inclination/position angle from the literature. The fixed model parameters are listed in Supplementary Table 2.

Outside rsnow, the CO line emissions arise from a warm layer above the mid-plane dust emission and thus the temperature of the dust continuum and the two CO isotopologues are significantly different. As previous studies have shown that the (sub)mm continuum emission is optically thin[38], we ignore dust emission in the outer disk. We use a constant temperature of 21 K for the CO isotopologues emission, as previously measured from CO isotopologue line emissions in the outer disk region of the TW Hya disk[39]. For the line optical depth, we use $\tau_{\rm C^{18}O}(r,v)=\delta\times\tau_{0}\times(r/r_{\rm snow})^{q_{2}}\times\left(1+Ae^{-\frac{(r-r_{\rm ring})^{2}}{2\sigma_{\rm ring}^{2}}}\right)\times e^{-r/100AU}\times\phi_{v}$, where a Gaussian ring is added to model the ring structure around 60 AU seen in the C18O (3-2) line emission[40, 39], and the exponential component ($e^{-r/100AU}$) is added to model the apparently sharp truncation of the CO emission beyond 100 AU.

For the continuum emission, we employ the radial profile of surface brightness derived from the high spatial resolution dust continuum observations by Andrews et al. (2016)[15], which were taken at a similar frequency (345 GHz) and with a much higher spatial resolution ($\sim$1 AU). We note that the spatial resolution of our constraints on the dust optical depth profile is 5-10 AU, because it is limited by the resolution of the temperature distribution measured from the C18O line emission (see Supplementary Figure 4).

To compare models with observations, we carry out Fourier transforms on the model image cubes and retrieve the complex visibilities at the ($u,v$)-plane locations sampled by the ALMA observations using the Python package vis_sample111Available at https://github.com/AstroChem/vis_sample. The modeled visibilities can then be directly compared with observed visibilities from ALMA observations. To sample the posterior probability distribution for each parameter, we use an MCMC approach. In practice, the sampling exploration is carried out by the affine-invariant routine EMCEE[41]. The prior distribution of parameters is assumed to be uniform across a wide parameter space. Our chains consist of 100 walkers and 1000 steps. The best-fit parameters and their uncertainties are listed in Table 1 and direct comparison of the observed C18O and 13C18O $J$=3-2 channel maps with synthetic data from the best model are plotted in Supplementary Figure 8 and 9.

0.3 CO snowline or jump in gas surface density.

To investigate whether the jump observed in the CO column density can be caused by a jump in the gas column density rather than the CO mid-plane snowline, we include two additional models of the gas mass distribution (Model 2 and 3). The two additional models are based on the assumption the gas-to-dust mass ratio and the dust opacity are constant throughout the disk. Model 2 is the gas surface density distribution derived by van boekel et al. (2016)[42] that fits the scattered light images of the TW Hya disk at 0.63, 0.79, 1.24 and 1.62 $\mu$m (expected to trace micron-sized dust grains that are coupled with the gas). Model 3 is based on the observed surface brightness distribution of the 345 GHz continuum emission of the TW Hya disk. Using the mid-plane temperature profile we measured from the C18O (2-1) line emission, we derive a surface density of the mm-sized dust particles and scale it by 100 to represent the gas mass distribution. The two models as well as our best-fitting model for the CO lines (Model 1) are plotted in the Supplementary Figure 10. For the three models, we show the radial surface brightness profiles of the C18O  (2-1) moment zero map and the line emissions at the 3.22 km/s velocity channel. We find that neither models based on dust observations (Model 2 and 3) could reproduce the observed C18O emissions. Therefore the CO column density jump cannot be explained by a jump in the gas surface density profile and thus chemical effect must be included. Considering that our best-fitting model requires an order of magnitude decreases in the CO column density and the temperature measured is consistent with the condensation temperature of CO onto water ice surface, the jump in CO gas column density is most likely caused by a mid-plane CO snowline.

0.4 Comparison with previous models.

The retrieved physical properties at the mid-plane show some differences compared to previous model predictions. We measured a mid-plane CO snowline at 20.5$\pm 1.3$ AU from the spatial resolved C18O and 13C18O (3-2) images, which is closer to the central star than the previous estimation of 30 AU that derived from the N2H*+* (4-3) line emission as an indirect tracer[12]. This discrepancy is not surprising as Van’t hoff et al. (2016) showed that N2H*+* line emissions provide only an upper limit to the mid-plane CO snowline radius[43]. Moreover, our CO snowline location is consistent with the prediction of 19 AU from the Van’t hoff et al. (2016) chemical models of the TW Hya disk. The mid-plane temperature at the CO snowline is 27${}^{+4}_{-3}$ K, consistent with the expected freeze-out temperature when CO condenses onto a surface formed by a mixture of water and CO2 ice[44].

We find that the dust emission has $\tau\sim$0.5 inside 20 AU, about a factor of 2-4 less than previous global models for the TW Hya disk[38]. The previous models, however, were based on continuum observations that did not resolve the inner 20 AU (beam size $\sim$20-25 AU) and a temperature profile derived from spectral energy distribution rather than direct measurements, which makes the estimation for the inner 20 AU zone highly uncertain. Moreover, if the dust emission is already optically thick, its brightness temperature at the CO snowline would be equal to the CO condensation temperature. However, Andrews et al. (2016)[15] showed that the brightness temperature is only 12 K at 20 AU, much cooler than realistic condensation temperature (21-46 K) in a high-pressure environment such as the disk mid-plane. Therefore, the dust emission inside 20 AU is possibly has smaller optical depth than previously assumed.

We measure a 12C/13C ratio of 40${}^{+9}_{-6}$ based on the optical depth ratio of the C18O/13C18O (3-2) line emissions. This ratio is smaller than the average ratio in local ISM of 69$\pm$6[13] but consistent with that expected in dense photo-dissociation regions[45]. The ratio is mostly constrained by the weak 13C18O emission on the west side of the disk and future higher sensitivity observations are needed to confirm this ratio. To test the sensitivity of parameters on the 12C/13C ratio, we also run models with a fixed 12C/13C ratio of 69. We find that the best-fitting parameters in these models are consistent with the results in Table 1 within uncertainties, except that the values of $\tau_{0}$ and $\delta$ increase by a factor of $\sim$2.

0.5 Gas column density.

With the radial profiles of temperature and line+dust optical depth from our parameterized model, we can derive the radial distribution of 13C18O column density by assuming that CO is in local thermal equilibrium (Supplementary Figure 7). This assumption is reasonable since the gas number density at the mid-plane is expected to be much higher than the critical density of the CO $J$=3-2 line. To further convert the 13C18O column density distribution to a total gas mass distribution, we use the spatially unresolved HD (1-0) line flux to estimate the 13C18O abundance[9]. The HD (1-0) line has been used as a total mass tracer of the TW Hya disk since the D/H ratio is well calibrated as 1.5$\times$10*-5* in objects within $\sim$100 pc of the Sun[46]. With an upper state energy of 128.5 K, the HD (1-0) line emission is dominantly from the warm region of the disk (T$>$20 K), similar to the warm region where CO remains in the gas phase. Therefore, our 13C18O gas column density distribution makes a good approximation to the mass distribution of HD that contributes to the HD (1-0) emission, given $x_{\rm{}^{13}C^{18}O}$ is a constant in the warm region throughout the disk. The HD (1-0) line flux is (6.3$\pm$0.7)$\times$10*-18* W m*-2* at 112 $\mu$m[9]. Using our 13C18O column density and temperature distribution, we find that 90% of the HD (1-0) line flux comes from the inner 20 AU of the disk, a region for which the dust continuum emission at 112 $\mu$m is optically thick ($\tau_{17{\rm AU}}\sim 3.8$ for mm-sized dust grains), and therefore the line flux arises only from one side of the disk. Following the approach in Bergin et al. (2013)[9], we assume the HD (1-0) line is in LTE and scale the 13C18O column density to match the measured HD line flux. Propagating uncertainties in our best-fitting parameters and in the HD (1-0) line flux, we derive an average 13C18O abundance of 1.7${}^{+1.3}_{-0.8}\times$10*-10* throughout the warm region of the TW Hya disk (T$>$20 K). This corresponds to a gas surface density profile of $\Sigma_{g}$ = 13${}^{+8}_{-5}\times$(R/20.5 AU)${}^{-0.87^{+0.38}_{-0.26}}$ g cm*-2* between 5-20.5 AU.

0.6 The gas-to-dust mass ratio and its uncertainty sources.

Combining the uncertainties in our gas and dust surface density distribution, we derive a total gas-to- dust ratio of 140${}^{+43}_{-31}$ for the 5-20.5 AU region. The main uncertainty for the gas mass estimation is from the gas temperature because the HD (1-0) line emission has a strong temperature dependence[9]. For the dust mass, we have not included the uncertainty of dust opacity due to dust growth. We use a $\kappa$ consistent with a dust population with an amax of 1 mm, but the $\kappa$ would decrease by a factor of two for an amax of 1 cm, and by a factor of 8 for 10 cm. For this reason, we treat our estimation as a gas-to-mm-sized dust mass ratio.

0.7 Sanity check.

The critical assumption in the model above is that the 13C18O, C18O, and dust continuum emission within the CO snowline can be characterized with the same temperature at each radius. Here we perform a sanity check to see if the assumption is valid in our best-fit model. The vertical distribution of gas and dust is characterized by their scale heights. Assuming a hydrostatic disk, we calculate the gas scale height using hg = $c_{s}/\Omega_{K}$, where $c_{s}$ is the disk sound speed at the mid-plane and $\Omega_{K}$ is the Keplerian angular velocity. The scale height of dust depends on turbulence level and the strength of the coupling between dust particles and turbulent gas[47]. Assuming a medium $\alpha$ value of 10*-3*, we find that the effective scale height of dust grains responsible for the 345 GHz continuum emission is $\sim$40% of that of the gas in the 5-21 AU for our best-fit model. As shown in Supplementary Figure 7, the 13C18O, C18O $J$=3-2 and dust optical depth is nearly flat inside 21 AU and can be characterized with constant values for each species. Using a $h_{d}/h_{g}$ of 0.4 and three typical optical depths for the 13C18O, C18O and continuum emission, we show that the optically thin 13C18O and dust emission are from both sides of the disk in vertical regions within two gas scale heights (Figure 2). Although the optically thick C18O line emission is mainly from one side of the disk, 80% of its flux are from the region within two scale heights. Therefore the best-fit model is self-consistent in the treatment of the temperature.

In more extreme cases that $\alpha$ approaches 10*-4*, the mm-sized grains become significantly more settled than the gas[29]. To check the temperature difference between gas (micron-sized dust) and mm-sized dust grains, we run radiative transfer models for the TW Hya disk with two dust populations. Both populations have a power-law size distribution n(a)$\propto$a*-3.5* between amin = 0.005 ${\rm\mu}$m and amax. The first population takes 10% of the total dust mass, has an amax=1 ${\rm\mu}$m, and thus is vertically couple with the gas. The second population accounts for 90% of the dust mass, has an amax=1 mm, and has a scale height of 20% of the gas scale height. We find that although the large grain population is cooler than the small grain population in the surface layer and outer disk region, the two populations share very similar temperature within one gas scale height between 5-20 AU (Supplementary Figure 11). In summary, it is reasonable to assume that 13C18O, C18O, and dust continuum emission within 20 AU can be characterized with the same temperature.

0.8 Data Availability Statement

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Supplementary Materials