<figure><img src="image_1.png" /></figure>

# Temperature, Mass, and Turbulence: A Spatially Resolved Multiband Non-LTE Analysis of CS in TW Hya

Richard Teague<fnref n="1" />,<fnref n="2" />, Thomas Henning<fnref n="2" />, Stéphane Guilloteau<fnref n="3" />, Edwin A. Bergin<fnref n="4" />, Dmitry Semenov<fnref n="2" />, Anne Dutrey<fnref n="3" />, Mario Flock<fnref n="2" />, Uma Gorti<fnref n="5" />, and Tilman Birnstiel<fnref n="6" />

Received 2018 June 7; revised 2018 July 31; accepted 2018 August 2; published 2018 September 10

## Abstract

Observations of multiple rotational transitions from a single molecule allow for unparalleled constraints on the physical conditions of the emitting region. We present an analysis of CS in TW Hya using the *J* = 7–6, 5–4 and 3–2 transitions imaged at ~0.5′′5 spatial resolution, resulting in a temperature and column density profile of the CS emission region extending out to 230 au, far beyond previous measurements. In addition, the 15 kHz resolution of the observations and the ability to directly estimate the temperature of the CS emitting gas, allow for one of the most sensitive searches for turbulent broadening in a disk to date. Limits of $v_{\text{turb}} \lesssim 0.1c_s$ can be placed across the entire radius of the disk. We are able to place strict limits of the local H₂ density due to the collisional excitations of the observed transitions. From these we find that a minimum disk mass of $3 \times 10^{-4} M_{\text{Sun}}$ is required to be consistent with the CS excitation conditions and can uniquely constrain the gas surface density profile in the outer disk.

**Key words:** astrochemistry – ISM: molecules – protoplanetary disks – techniques: interferometric

## 1. Introduction

To understand the planet formation process we must first understand the environmental conditions of planetary birth-places (Mordasini et al. 2012). Thanks to the unparalleled sensitivity and resolution provided by the Atacama Large (sub-)Millimeter Array (ALMA), we are routinely resolving substructures indicative of in situ planet formation and on-going physical processing (ALMA Partnership et al. 2015; Andrews et al. 2016; Pérez et al. 2016; Fedele et al. 2017; Dipierro et al. 2018). Similar features have been observed in high-contrast AO near-infrared imaging tracing the small grain population well coupled to the gas (e.g., van Boekel et al. 2017; Pohl et al. 2017; Hendler et al. 2018).

In addition, excess UV emission is interpreted as accretion onto the central star from the disk, with more massive disks accreting at a higher rate (Fang & White 2004; Manara et al. 2016). These observations point toward active disks that are able to efficiently redistribute material and angular momentum.

Despite the observational evidence for the redistribution of angular momentum, the physical mechanisms that enable this remain elusive. Two scenarios are likely. First, a turbulent viscosity would be sufficient to transport angular momentum outward. This is the assumption in the frequently implemented “α-viscosity” disk model of Shakura & Sunyaev (1973). Although agnostic about the main driver of the turbulence, this model links the turbulent motions to the viscosity of the disk. Alternatively, angular momentum can be removed through winds (Bai 2017), evidence of which has been observed in several young sources, though it is lacking more evolved counterparts.

The magnetorotational instability has been a leading contender as the source for turbulence (Balbus & Hawley 1998;

Fromang & Nelson 2006; Simon et al. 2013, 2015; Bai 2015; Flock et al. 2015, 2017). However, estimates of the local ionization rate close to the disk midplane have suggested that there would be insufficient coupling between the rotating gas and the magnetic field (Cleeves et al. 2015). Additional instabilities have been shown to generate turbulence without the need for ionization such as the vertical shear instability (Nelson et al. 2013; Lin & Youdin 2015), gravitational instabilities (Gammie 2001; Forgan et al. 2012), baroclinic instabilities (Klahr & Bodenheimer 2003; Lyra & Klahr 2011), and the zombie vortex instability (Marcus et al. 2015; Lesur & Latter 2016). Distinguishing between these mechanisms requires a direct comparison of the distributions of nonthermal motions observed in a disk and the predicted distribution from simulations (for example, Forgan et al. 2012; Flock et al. 2015; Simon et al. 2015).

There have been several attempts to detect nonthermal motions in disks through the additional broadening in line emission in the disks of TW Hya and HD 163296 (Hughes et al. 2011; Guilloteau et al. 2012; Flaherty et al. 2015, 2017, 2018; Teague et al. 2016). Although a promising avenue of exploration, this approach is hugely sensitive to the temperature assumed, as the Doppler broadening of the lines does not distinguish between the sources of the motions, either thermal or nonthermal. One must make assumptions about the physical structure of the disk in order to break these degeneracies (Flaherty et al. 2015, 2017; Simon et al. 2015). Teague et al. (2016) attempted to minimize the assumptions made about the disk structure, inferring physical properties directly from the observed spectra and allowing the temperature and turbulent structure to vary throughout the disk. Without the leverage of an assumed model, the constraints on $v_{\text{turb}}$ were larger than other attempts, finding $v_{\text{turb}} \lesssim 0.3 c_s$ across the radius of the

<footnote n="1">¹Department of Astronomy, University of Michigan, 1085 S. University Avenue, Ann Arbor, MI 48109, USA; rteague@umich.edu</footnote>

<footnote n="2">²Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="3">³Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, Alle Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="4">⁴Department of Astronomy, University of Michigan, 311 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA</footnote>

<footnote n="5">⁵SETI Institute/NASA Ames Research Center, Mail Stop 245-3, Moffett Field, CA 94035-1000, USA</footnote>

<footnote n="6">⁶University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 Munich, Germany</footnote>
<table><caption>**Table 1**
List of Observations</caption><thead><tr><th rowspan="2">CS Transition</th><th rowspan="2">Frequency (GHz)</th><th rowspan="2">E<sub>u</sub> (K)</th><th rowspan="2">ΔV<sub>chan</sub> (m s<sup>-1</sup>)</th><th rowspan="2">Robust</th><th rowspan="2">UV Taper (kλ)</th><th colspan="2">Restoring Beam</th><th rowspan="2">Image rms<sup>a</sup> (K)</th><th rowspan="2">Peak T<sub>B</sub> (K)</th></tr><tr><th>(" × ")</th><th>(°)</th></tr></thead><tbody><tr><td>J = 3–2</td><td>146.96904</td><td>14.11</td><td>35</td><td>2.0</td><td>320</td><td>0.57 × 0.50</td><td>91.3</td><td>0.56</td><td>10.7</td></tr><tr><td>J = 5–4</td><td>244.93556</td><td>35.27</td><td>19</td><td>0.5</td><td>...</td><td>0.59 × 0.47</td><td>91.4</td><td>0.46</td><td>11.6</td></tr><tr><td>J = 7–6</td><td>342.88285</td><td>65.83</td><td>14</td><td>2.0</td><td>230</td><td>0.57 × 0.51</td><td>74.9</td><td>0.39</td><td>10.3</td></tr></tbody></table>

**Note.**

<sup>a</sup> The resulting noise for the azimuthally averaged spectra is reduced by a factor of $\sqrt{N}$, where N is the number of beams averaged over.

disk. These constraints were limited by assumptions made about the thermalization of the energy levels, in particular, CN emission was shown to be in non-LTE across the outer disk, while for the single CS transition, as the line was optically thin, the degeneracy between column density and temperature could not be broken without assuming an underlying physical structure.

However, high-resolution observations of the studied disks show substructures traced in millimeter continuum, molecular line emission and scattered light (Andrews et al. 2016; Flaherty et al. 2017; Monnier et al. 2017; Teague et al. 2017; van Boekel et al. 2017). Such perturbations from a “smooth” disk model could be sufficient to mask any signal from nonthermal motions, which require a precise measure of the temperature to a few kelvin.

In this paper, we present new ALMA observations of CS J = 7–6 and J = 3–2 transitions in TW Hya, the nearest protoplanetary disk at $d = 60.1$ pc (Bailer-Jones et al. 2018) with a near face-on inclination of $i \approx 6^\circ$. Combined with the previously published J = 5–4 observations (Teague et al. 2016), we are able to fit for the excitation conditions of the molecule, namely the temperature, density, column density, and nonthermal broadening component.

## 2. Observations

### 2.1. Data Reduction

The new observations were part of the ALMA project 2016.1.00440.S, targeting the (7–6) and (3–2) rotational transitions of CS at 342.882850 GHz and 146.969029 GHz in Bands 7 and 4, respectively. Band 6 data come from project 2013.1.00387.S, originally published in Teague et al. (2016).

Band 4 data were taken in 2016 October 22, 25, and 27 with 38, 39, and 40 antennae, respectively, with baselines spanning 18.58 m–1.40 km and a total on-source time of 138.73 minutes. Band 7 data were taken on the 2016 December 2 utilizing 45 antennas with baselines spanning 15.10 m to 0.70 km and a total on-source time of 47.75 minutes. The quasars J1107–449 and J1037–2934 were used as amplitude and phase calibrators for both bands. For each band, the correlator was set up to have a spectral window covering the CS transition with a channel spacing of 15.259 kHz.

The data were initially calibrated with the provided scripts. Phase self-calibration was performed using CASA v4.7. Gain tables were generated on collapsed continuum images, then applied to the line spectral windows. Imaging was performed in CASA v5.2. Continuum was subtracted from all line spectral windows using the uvcontsub task. As continuum emission is observed to extend out to $\approx 70$ au (Andrews et al. 2016), this subtraction will not affect the line emission outside this region. Within 70 au, however, regions where the continuum is

optically thick may suffer from small absorption effects (Boehler et al. 2017); however, this will be limited to the very inner regions, which are not considered in this work.

All cubes were CLEANed using a Keplerian mask with parameters consistent with previous observations of TW Hya. To mitigate any differences between observations due to the different beam sizes, we applied a different weighting scheme to each transition. With the least extended baselines, the band 6 data were imaged using a robust weighting using a Briggs parameter of 0.5. To attain a comparable beam size with the band 4 and 7 data, we used natural weighting and included a Gaussian taper to the extended baselines. Each image was centered on the peak of the continuum using the fixvis command, with the centers being consistent within 50 mas in both right ascension and declination. Table 1 summarizes the observations and resulting images.

The ALMA Technical Handbook quotes a flux calibration uncertainty of 5% for Band 4 and 10% for Bands 6 and 7. Comparison of the integrated continuum fluxes for our data with other published values suggests that there are no significant deviations, as discussed in Appendix A. In the following, we assume a 10% flux calibration uncertainty for all three bands.

### 2.2. Observational Results

Zeroth and first moment maps were generated using both the Keplerian mask used for cleaning and clipping values below $2\sigma$. These are shown in Figure 1. All three lines display a similar emission pattern with an off-center peak at $\sim 1''$ (60 au) and extending to $\sim 3''$ (180 au). There is no clear azimuthal asymmetry within the noise. In addition, all lines exhibit the characteristic pattern of Keplerian rotation.

To derive geometrical properties for the disk, we fit a Keplerian rotation pattern, $v_{Kep} = \sqrt{GM_*/r}$ to the first moment maps. The disk center, $(x_0, y_0)$, inclination $i$, position angle, PA and systemic velocity, $v_{LSR}$, were allowed to vary. Following Walsh et al. (2017), we convolve the rotation pattern before calculation of the likelihood. The stellar mass was fixed at $0.65 M_{\text{Sun}}$ (Qi et al. 2004) because the disk inclination is too low to break the degeneracy between $M_*$ and $\sin(i)$ in the velocity pattern.

Results for each line were consistent, yielding an inclination $i = 6^\circ.8 \pm 0^\circ.1$, PA = $151^\circ \pm 1^\circ$, and $v_{LSR} = 2842 \pm 2 \text{ ms}^{-1}$. The uncertainties on these values are dominated by the scatter between the three observations and the statistical uncertainties are an order of magnitude smaller. The derived $M_* \cdot \sin i$ value is consistent with the $0.88 M_{\text{Sun}}$ and $i = 5^\circ$ used to model high-resolution CO emission (Huang et al. 2018). The use of these values would not affect the results in the following analysis.
<figure><img src="image_2.png" /><figcaption>Figure 1. Moment maps of the three CS emission lines: zeroth moment, top, and first moment, bottom. Beam sizes are shown by the hatched ellipse in the bottom left of each panel.</figcaption></figure>

Following Teague et al. (2016; also see Yen et al. 2016; Matrà et al. 2017), we deproject the data spatially and spectrally accounting for the rotation of the disk. Each pixel at deprojected coordinate (r, θ) was shifted by an amount

$$v_{\text{proj}}(r, \theta) = v_{\text{Kep}}(r) \cdot \sin(i) \cdot \cos(\theta), \quad (1)$$

to a common velocity. The data were then binned into annuli with a width of 4 au and then averaged to increase the signal-to-noise ratio. Although this binning is much below the beam size, broader annuli would sample lines arising from significantly different temperatures and thus possessing different line widths. Although in the following text we treat each annulus as independent, in practice, they will still be correlated over the size of the beam.

A Gaussian line profile, characterized by a line width ΔV and the peak brightness temperature T<sub>B</sub>, was fit to each azimuthally averaged spectrum (as shown in Appendix B). A

Gaussian profile is a reasonable assumption for such lines that are believed to dominated by thermal broadening. As the CS emission is expected to arise in a narrow layer (Dutrey et al. 2017), we do not expect significant deviations in the rotation velocity that would lead to non-Gaussian profiles. Furthermore, the face-on orientation of TW Hya means that the near and far sides of the disk align along the line of sight minimizing the contamination found in more inclined disks (for example, HD 163296, Rosenfeld et al. 2013), suggesting that a Gaussian profile is a reasonable assumption.

The radial profiles of these parameters are shown in Figure 2. The broad ring of emission centered at 60 au is clearly seen in the T<sub>B</sub> profile, while the knee at 90 au is seen in all three transitions. Teague et al. (2017) argued that this feature was the result of a surface density perturbation. An additional feature in the outer disk, either a drop in T<sub>B</sub> at ≈170 au or an increase at ≈190 au, is also seen, potentially an outer desorption front from interstellar radiation.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Radial profiles of the brightness temperature, $T_B$, and the line width $\Delta V$ derived by fitting a Gaussian profile to the deprojected spectra. All uncertainties are $3\sigma$ statistical uncertainties from the fit.</figcaption></figure>

Brightness temperatures peaking at $T_B \approx 9$ K suggests that the lines are optically thin. Assuming that the line width in the outer disk is dominated by thermal broadening, we expect temperature of $T_{kin} \gtrsim 25$ K, suggestive of optical depths of $\tau \lesssim 0.4$. This temperature is consistent with the freeze-out temperature of CS, indicating that the CS emission traces the CS freeze-out surface, comparable to that observed for the optically thin CO isotopologues (Schwarz et al. 2016).

### 3. LTE Analysis

In this section, we assume that the three observed transitions are in local thermodynamic equilibrium (LTE) so that $T_{ex} = T_{kin}$. The excitation temperature will provide a lower limit to gas temperature from the three lines. Loomis et al. (2018) recently used the same approach to infer the temperature and column density of methyl cyanide, CH$_3$CN, in TW Hya.

Following Goldsmith & Langer (1999), in the optically thin limit the integrated flux can be related to the level population of the upper energy state,

$$ \frac{N_{\text{u}}}{g_{\text{u}}} = \frac{1}{g_{\text{u}}} \frac{8\pi k \nu^2}{A_{\text{ul}} h c^3} \cdot \int_{T_B} dV = \frac{\gamma_{\text{u}} W}{g_{\text{u}}}, \quad (2) $$

where $W = \int_{T_B} dV$ is the value described by the zeroth moment maps in Figure 1, transformed to units of kelvin using Planck's law, $N_u$ is the column density of the upper energy level with degeneracy $g_u$, $A_{ul}$ is the Einstein-A coefficient for spontaneous emission, and $\nu$ is the frequency of the emission. Under the assumption of LTE, we can relate

$$ \ln \left( \frac{\gamma_u W}{g_u} \right) = \ln N - \ln Q(T) - \ln C_{\tau} - \frac{E_u}{k T_{\text{ex}}}, \quad (3) $$

where $N$ is the total column density of the molecule, $Q$ is the partition function, approximated for a linear rotator as,

$$ Q(T) = \frac{kT}{hB_0} + \frac{1}{3}, \quad (4) $$

and $C_\tau$ is the optical depth correction factor, $C_\tau = \tau/(1 - \exp(-\tau))$ for a square line profile. In the case of a Gaussian profile, this correction is less severe. The optical depth at the line center is given by

$$ \tau = \frac{N_u A_{ul} c^3}{8\pi \Delta V_{\text{FWHM}} \nu^3} \left[ \exp\left(\frac{h\nu}{kT_{\text{ex}}}\right) - 1 \right], \quad (5) $$

where $\Delta V_{\text{FWHM}}$ is the full width at half maximum of the line, which can be calculated assuming only thermal broadening, a reasonable assumption given the low levels of turbulence previously reported.

Finally, relating $N_u$ to $N$ through

$$ N_u = \frac{g_u N}{Q(T)} \exp\left(\frac{-E_u}{kT_{\text{ex}}}\right), \quad (6) $$

allows us to fit for {$T_{ex}, N(\text{CS})$}, self-consistently accounting for possible optical depth effects. As this method considers only the integrated flux, flux calibration uncertainties can be easily considered. We include a systematic flux calibration uncertainty of 10% in addition to the statistical uncertainty.

As this approach uses only the integrated flux value there are a range of {$T_{ex}, N(\text{CS})$} that are consistent with the data but result in highly optically thick lines. We can rule these scenarios out given the $T_B$ profiles shown in Figure 2: optically thick lines result in $T_B = T_{ex}$, recovering a gas temperature far below the freeze-out temperature of volatile species. To take this into account, we include a prior that $\tau < 1$ for all lines.

We use emcee (Foreman-Mackey et al. 2013) to calculate posterior distributions of $T_{ex}$ and $\log_{10} N(\text{CS})$. We assumed uninformative priors,

$$ \begin{aligned} T_{\text{ex}} (\text{K}) &= \mathcal{U}(10, 150) \\ \log_{10}(N(\text{CS}) (\text{cm}^{-2})) &= \mathcal{U}(9, 14) \\ \tau (T_{\text{ex}}, N(\text{CS})) &= \mathcal{U}(0, 1), \end{aligned} \quad (7) $$

forcing the models to be optically thin. We used 256 walkers that took 200 steps to burn-in before taking an additional 100 to sample the posterior distribution. The quoted uncertainties are the 16th and 84th percentiles of the posterior distribution, which for a Gaussian distribution represents $1\sigma$ uncertainty.

The results are shown in Figure 3. The left two panels show the 2D maps of $T_{ex}$ and $\log_{10} N(\text{CS})$ where the fitting was applied on a pixel-by-pixel basis. The two panels on the right show the results when applied to the azimuthally averaged spectra. The blue dots include the correction for optical depth while the gray dots do not. It is only within the inner regions where $\tau$ rises where the correction for the optical depth is appreciable.

Within the inner 30 au a constant temperature of $T_{ex} \approx 40$ K is found, slowly dropping to 20 K at 200 au. The column density shows two distinct knees at 90 au and 160 au, the former argued for in Teague et al. (2017) as due to a surface density perturbation. A slight rise in temperature is observed at 150 au.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** LTE analysis for the three CS transitions. Left and center panels show $T_{ex}$ and $\log_{10} N(\text{CS})$ applied on a pixel-by-pixel basis. The beam in the bottom left shows the average beam of the three transitions. The right two panels show the LTE approach applied to the azimuthally averaged spectra. Blue dots show the results including the optical depth correction, $C_\tau$, and red show without. The error bars show the 16th to 84th percentiles of the posterior distributions.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

The 2D distributions show no significant azimuthal structure in either parameter with deviations being consistent with the uncertainty and show a comparable radial profile to the azimuthally averaged spectra. The 2D maps are only able to achieve a reasonable fit within ≈2" because of the noise, demonstrating the strength of the azimuthally averaged approach in tracing the outer disk.

## 4. Non-LTE Modeling

By collapsing our data down to a single integrated flux value as in the previous section, we lose a tremendous amount of information and modeling the entire spectrum allows us to consider a more complex model. In this section, we fit the spectra using slab models, allowing us to explore the impact of non-LTE effects and to potentially constrain the local gas density.

We assume that the line profile is well described by an isothermal slab model with

$$T_B(\nu) = (J_\nu(T_{ex}) - J_\nu(T_{bg})) \cdot (1 - e^{-\tau_v}) + J_\nu(T_{bg}), \quad (8)$$

where

$$J_\nu(T) = \left( \frac{h\nu/k}{\exp(h\nu/kT) - 1} \right), \quad (9)$$

$T_{bg}$ = 2.73 K is the background temperature and

$$\tau_v = \tau_0 \exp\left(-\frac{(v - v_0)^2}{\Delta V^2}\right), \quad (10)$$

is the optical depth (Rohlfss & Wilson 1996). Such a profile accounts for the saturation of the line core at moderate optical depths where the line profile can deviate significantly from a Gaussian profile (see, for example, Teague et al. 2016). The line width is the quadrature sum of thermal broadening and nonthermal broadening components,

$$\Delta V = \sqrt{\frac{2kT_{\text{kin}}}{\mu m_p} + v_{\text{turb}}^2}, \quad (11)$$

where $\mu = 44$ is the molecular weight of CS and $v_{\text{turb}}$ is the line-of-sight velocity dispersion. Nonthermal broadening is parameterized as a fraction of the local sound speed, $M = v_{\text{turb}}/c_s$, where $c_s$ is the local sound speed. This allows for a rough comparison with the frequently used $\alpha$ viscosity parameter via the relation $\alpha \sim M^2$. However, as this relation is dependent on the form of the viscosity (Cuzzi et al. 2001), we limit our discussion to only the Mach number.

From a gas kinetic temperature $T_{\text{kin}}$, the local gas volume density $n(\text{H}_2)$ and the column density of CS, $N(\text{CS})$, we can calculate the excitation temperatures for each observed level (where $T_{\text{ex}} \le T_{\text{kin}}$) and optical depth at the line center, $\tau_0$. We use the collision rates from Denis-Alpizar et al. (2018) and assume a thermal $\text{H}_2$ ortho-para ratio,

$$\frac{n(\text{ortho} - \text{H}_2)}{n(\text{para} - \text{H}_2)} = 9 \times \exp\left(-\frac{170.5 \text{ K}}{T_{\text{kin}}}\right) \quad (12)$$

from Flower & Watt (1984, 1985). For a gas of 30 K, this is ~0.03. These values can then be used to model a full line profile through Equations (8) and (10).

In practice, the excitation calculations are performed with the 0D code RADEX (van der Tak et al. 2007), assuming a slab geometry, appropriate for the face-on viewing orientation of TW Hya. To achieve the speed necessary for MCMC sampling, a large grid in {$T_{\text{kin}}$, $n(\text{H}_2)$, $N(\text{CS})$, $\Delta V$} space was run using pyradex,<fnref n="7" /> and linearly interpolated between points. Both $n(\text{H}_2)$ and $N(\text{CS})$ were sampled logarithmically while $T_{\text{kin}}$ and $\Delta V$ were linearly sampled, with 40 samples along each axis. The resulting grid was checked to confirm that $\tau$ and $T_{\text{ex}}$ were smoothly varying over the grid ranges and that linear interpolation was appropriate.

As we are working in the image plane, before comparing the synthetic spectra with data, we must make corrections for the sampling of the ALMA correlator and the imaging process. Flaherty et al. (2018) discuss the difference in image-plane and uv-plane fitting for interferometric data, concluding that for observations with well sampled uv-planes, such as the data

<footnote n="7">⁷ https://github.com/keflavich/pyradex</footnote>
presented here, there is little difference in the results. Nonetheless, caution must be exercised as spatial correlations could lead to underestimation of the uncertainties.

Spectra are generated at a sampling rate of 150 Hz, a sampling rate of 100 times the observations, before being sampled down to 15 kHz and Hanning smoothed by a kernel with a width 15 kHz. This step is essential, as with such narrow line widths this can result in underestimating the intrinsic peak brightness temperature by ~20% and overestimating the width by ~10%. Figure 13 from Rosenfeld et al. (2013) demonstrates how this process affects the emission morphology.

In addition to modeling the emission, we attempt to account for possible spectral correlations in the noise by modeling the noise component using Gaussian Processes, implemented with the Python package celerite (Foreman-Mackey et al. 2017). A Gaussian Process is a probabilistic nonparametric approach to modeling smoothly varying functions and readily allows for the inclusion of covariances between data points. The noise was modeled with an approximate Matern-3/2 kernel, an approximation of a Gaussian to mimic the Hanning smoothing applied to the data in the correlator, which is specified by the amplitude of the noise and the correlation length in velocity, $\sigma_{rms}$ and $\ell$, respectively. The kernel describes how points are correlated over a given dimension, in this case the spectra axis. Figure 1 of Czekala et al. (2017) provides and example of such kernels. A simple harmonic oscillator kernel was also tested; however, no significant difference was found between kernels. These parameters were ultimately considered nuisance parameters, used only to consider how different noise models affected the results, and marginalized over in the calculation of the excitation condition posterior distributions.

Due to the finite resolution of the data, the synthesized beam will smear out the emission spatially. As discussed in Teague et al. (2016), this will lead to a broadening of the lines in the inner region of the disk, while at larger radii, these effects are negligible. We do not attempt to correct for such beam smear effects. The effects are two-dimensional in nature and attempting to model them in one dimension introduces significant uncertainty to the modeling procedure.

In total, for each radial position the three lines can be fully specified by 13 free parameters: the local physical conditions, {$T_{kin}$, $n(H_2)$, $N(CS)$, $M$}; the center of each line, {$v_0$}; and the noise model for each observation {$\sigma$, $\ell^2$}. All parameters were given uninformative (uniform) priors, ranging across values expected in a protoplanetary disk,

$$ 
\begin{aligned}
T_{\text{kin}} (\text{K}) &= \mathcal{U}(10, 150) \\
\log_{10}(n(\text{H}_2) (\text{cm}^{-3})) &= \mathcal{U}(4, 10) \\
\log_{10}(N(\text{CS}) (\text{cm}^{-2})) &= \mathcal{U}(9, 14) \\
\log_{10} M &= \mathcal{U}(-5, 0).
\end{aligned} 
\quad (13) 
$$

We additionally run a set of models where we fit for $M$ rather than $\log_{10} M$ for which we impose a prior of $M = \mathcal{U}(0, 1)$. Each fit consisted of 256 walkers, each taking 1500 steps for a burn-in period, then an additional 250 to sample the posterior distribution.

The results are shown in Figure 4, with blue points showing the fit to $\log_{10} M$ and gray points to the linear fit. Example of the covariances between parameters and their posterior distributions can be found in Appendix C. The error bars denote the 16th to 84th percentiles of the posterior distributions

<figure><img src="image_6.png" /><figcaption>Figure 4. Results of the non-LTE fitting. Blue values show the fits to $\log_{10} M$, while the red points show fits to $M$. The error bars show the 16th to 84th percentile range of the posterior distribution and the central point is the median. For panel (d), triangles show upper limits for the value, defined as when the 16th percentile of $M \le 10^{-2}$.</figcaption></figure>

around the median value. Both approaches yield comparable values, consistent with the LTE approach described in Section 3. This is not surprising as high $H_2$ densities are found resulting in thermalization of the transitions out to $\approx 190$ au.

We see a significant deviation between the models, one where $M$ is varied linearly and the other where it is varied logarithmically, in the inner 90 au, which is due to the artificial broadening of the line from the beam smear. In brief, the broadened lines require either a larger temperature or a large nonthermal broadening component. When fitting for $\log_{10} M$, changes in the temperature are preferred over changes in the nonthermal broadening due to the temperature only being considered linearly. Conversely, the fit for $M$ prefers solutions with larger nonthermal broadening without requiring an increase in the temperature. At these higher temperatures assumed in the logarithmic fit, the $J = 3-2$ line is considerably less emissive than the $J = 7-6$ and thus to maintain a higher $J = 3-2$ line flux, the volume density must be artificially decreased to reduce the importance of collisional excitation.

The noise models were able to converge, resulting in Gaussian distributions of their parameters which were not
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparing a power-law profile fit to the CS column density in gray to the derived values in blue. The bottom panel shows the relative residual between the two.</figcaption></figure>

correlated with any other parameter. We have marginalized over these distributions when analyzing the distributions of the parameters describing the disk physical conditions. See Appendix C for an example of the posterior distributions.

# 5. Discussion

## 5.1. Disk Physical Structure

Both LTE and non-LTE approaches paint the same picture: CS is present across the entire extent of the disk in a region that slowly cools from ≈40 K in the inner disk to ≲20 K at 200 au.

One interpretation for this is that the CS layer is bounded by the CS “snow surface” (Schwarz et al. 2016; Loomis et al. 2018). The binding energy of CS is 1900 K resulting in a desorption temperature of $T_{desorb} \approx 31$ K, although dependent on local gas pressures that are unknown, comparable to the temperatures observed for CS. However, a clear boundary between gaseous and ice forms of CS may not be present due to the chemical reprocessing expected on the grain surfaces. Observations of high inclination disks would be able to identify whether the CS emission has a sharp lower boundary.

The column density profiles show two distinct knee features at 90 au and 160 au, both seen in the $T_B$ profiles in Figure 2. Figure 5 compares the column density derived in Section 4 with a power-law profile fit (shown by the black line). The residuals, shown in the lower panel, show deviations of up to 20% in $N(\text{CS})$ at 120 and 160 au. Despite these deviations, models assuming a simple power-law column density, such as those in Teague et al. (2016), would be able to adequately model the true column density profile.

The dip at 90 au was previously argued by Teague et al. (2017) to be due to a significant perturbation in the gas surface density needed to account for a gap traced by the scattered light (van Boekel et al. 2017). While these results confirm that the emission feature is due to a change in column density rather than temperature, they are unable to distinguish between a local change in CS abundance and a total depletion of gas. Similar features have been observed in high-resolution $^{12}$CO observations (Huang et al. 2018).

At the outer edge of the disk, the density drops to a sufficiently low value that, unlike inward of ≈190 au, the

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: upper limits on $\Delta V$ compared to the observed values that span the gray shaded region. Bottom: $2\sigma$ upper limits on the nonthermal broadening. The shaded regions show the limits from Teague et al. (2016) and Flaherty et al. (2018).</figcaption></figure>

volume density of $H_2$ can be constrained. The apparent constraints inward of 90 au are, as discussed before, an artifact of the limited angular resolution. Observations of higher frequency lines with higher critical densities will allow for this method to be sensitive to the higher densities of the inner disk and allow us to extend these surface density constraints to smaller radii.

## 5.2. Turbulence

Determination of the nonthermal broadening requires an accurate measure of the local temperature in order to account for the thermal contribution. By constraining the temperature through multiple transitions minimizes assumptions about the thermal structure and provides the most accurate measure of the gas temperature to date. As the CS lines are optically thin, this temperature will be the contribution function-weighted average of the emitting column. With our derived temperature profile, we are therefore able to derive spatially resolved limits on the required nonthermal broadening to be consistent with the data.

For TW Hya, multiple studies have already been undertaken (Hughes et al. 2011; Teague et al. 2016; Flaherty et al. 2018), finding a range of nonthermal broadening values and upper limits, $M \lesssim 0.4$. As discussed in Flaherty et al. (2018), differences in these limits are primarily driven by the different assumptions about the underlying thermal structure and how this couples to the density structure. Teague et al. (2016) caution, however, that constraints of $M \lesssim 0.03$ require one to constrain the thermal structure to near kelvin-precision, a limit that is achieved with the data presented in this manuscript.

Under the assumption that the presented three CS lines arise from the same vertical layer in the disk, an assumption that requires observations of edge-on disks to test, we are able to remain agnostic about the thermal and physical structure of the disk. As shown in the bottom panel of Figure 6, we are able to place a $2\sigma$ upper limit out to 230 au. Both approaches (linear and logarithmic fits of $M$ in red and blue, respectively) yield
comparable results. The rise in limits inward of 100 au is due to the broadening arising due to the beam smearing, while outside 180 au the limits increase due to the lower SNR of the data. The fits yields limits consistent with the values found by Flaherty et al. (2018), $M \le 0.13$. These are a factor of a few lower than Teague et al. (2016) due to the warmer temperature derived ($T_{kin} = 28$ K at 100 au compared to 12 K as in Teague et al. 2016) as only a single CS transition was available.

Two dips are observed in the profiles at 110 au and 165 au, consistent with the dips in the column density. The resulting line widths, plotted in the top panel of Figure 6 show that they yield comparable widths to the observations, however, overproduce the line width outside 200 au, which is likely due to the lower SNR of the data. We leave the interpretation of these features for future work.

Deriving limits for $M$ using a parametric modeling approach, where physical properties and chemical abundances are described as analytical functions, as in Guilloteau et al. (2012), Flaherty et al. (2015, 2017, 2018), and Teague et al. (2016) requires a well constrained molecular distribution that, for CS, is not known. Furthermore, there is mounting evidence that spatially varying abundances of C and O within the disk can radially alter the local chemistry (Bergin et al. 2016), further limiting the accuracy of a simple analytical prescription.

### 5.3. Minimum Disk Mass

Many studies have attempted to measure the mass of the TW Hya disk through observations of both the millimeter continuum and gas emission lines. With differing assumptions, these have resulted in a range of masses spanning $5 \times 10^{-4} M_{Sun}$ to $6 \times 10^{-2} M_{Sun}$ (Kastner et al. 1997; Calvet et al. 2002; Thi et al. 2010; Gorti et al. 2011; Favre et al. 2013). Arguably the most accurate approach is to use hydrogen deuteride, HD, as this molecule should be tracing the $H_2$ gas most closely. Modeling the HD $J = 1-0$ transition, Bergin et al. (2013) concluded that the mass of the disk must be $M_{disk} > 0.05 M_{Sun}$. More recently, Trapman et al. (2017) used additional observations of the $J = 2-1$ transition to find a mass between $6 \times 10^{-3}$ and $9 \times 10^{-3} M_{Sun}$. This large range is due primarily to the sensitivity of the HD emission to the assumed thermal structure and differences in assumed cosmic D/H ratio. Models of HD emission show that it is almost entirely confined to the warm inner disk, $r < 100$ au, where the gas is warm enough to sufficiently excite the fundamental transition. Although this region accounts for almost all the disk gas mass, the HD flux is insensitive to the cold gas reservoir at smaller radii and is thus a minimum disk mass.

As we have limits on the required $H_2$ density as a function of radius in Section 4, we are able to place a limit on the minimum gas surface density and thus disk mass needed to recover the inferred excitation conditions. It is important to note that this technique does not require the assumption of a molecular abundance, such as those using HD, but rather constrains the $H_2$ gas directly through collisional excitation. It therefore provides an excellent comparison for techniques that aim to reproduce emission profiles and provides a unique constraint for surface densities at large radii in the disk.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Inferred minimum gas surface density and resulting minimum disk masses using the derived $n(H_2)$ values. The blue line shows the results when fitting for $\log_{10} M$ and the red when fitting for $M$. The gray shaded region shows where $n(H_2)$ has been constrained. The dotted and dotted-dashed lines show the surface density profiles used to model the HD emission in Bergin et al. (2013; using the surface density from Gorti et al. 2011) and Trapman et al. (2017). The dashed line shows the surface density used by van Boekel et al. (2017) to model the scattered light emission profile.</figcaption></figure>

To scale a midplane density to a column density, we assume a Gaussian vertical density structure,

$$ \rho_{\text{gas}} = \frac{\Sigma_{\text{gas}}}{\sqrt{2\pi} H_{\text{mid}}} \times \exp\left(-\frac{z^2}{H_{\text{mid}}^2}\right), \quad (14) $$

where we take the pressure scale height,

$$ H_{\text{mid}} = \sqrt{\frac{kT_{\text{mid}}r^3}{\mu m_H G M_{\star}}}, \quad (15) $$

which is dependent on the assumed midplane temperature, $T_{mid}$. Observations of the edge-on Flying Saucer have shown CS emission to arise from $\lesssim 1H_{mid}$ (Dutrey et al. 2017). The angular resolution of these observations does not allow for distinction between the case of two elevated, thin molecular layers at $\pm 1H_{mid}$ or a continuous distribution below $H_{mid}$. Measurements of the midplane temperature estimate this to be 5–7 K for the millimeter dust (Guilloteau et al. 2016) and $\approx 12$ K for the gas (Dutrey et al. 2017). As we find $T_{kin} = 20–35$ K, this suggests that CS is not tracing the midplane, but rather a slightly elevated region, so that we would overestimate the pressure scale height and underestimate the midplane density. In combination, these uncertainties should mitigate one another allowing for a first-order estimation of the minimum surface density.

For this estimate, we use both fits from Section 4 for the minimum $n(H_2)$ to infer a minimum $\Sigma_{gas}$, shown in Figure 7. Integrating these minimum surface densities, we find an average minimum disk mass of $3 \times 10^{-4} M_{Sun}$, fully consistent with the estimates from HD emission. Observations of transitions that thermalize at higher densities would extend the sensitivity of this approach such that it can distinguish between models predicting different disk masses.

The shaded region at $r > 190$ au highlights the region where $n(H_2)$ was measured rather than a lower limit constrained. In this region, we expect the plotted minimum surface densities to be close to the $\Sigma_{gas}$ value rather than just a minimum value; however, the accuracy will be limited by the assumptions made
above the vertical structure of the disk. Nonetheless, these profiles provide unique constraints on the gas surface density in the outer regions and are highly complimentary to studies using optically thin CO isotopologues that trace Σ<sub>gas</sub> within the CO snowline (Schwarz et al. 2016; Zhang et al. 2017).

In Figure 7, we additionally plot surface density profiles from Gorti et al. (2011), used in Bergin et al. (2013) to model the HD flux, the best-fit profile from Trapman et al. (2017), also used to model the HD flux and finally the profile from van Boekel et al. (2017) used to model scattered light emission. From our lower limit, we are able to rule out the model from Gorti et al. (2011) that contains insufficient material in the outer disk to recover the excitation conditions required by the CS transitions. The profile from Trapman et al. (2017) is broadly consistent with the lower limits; however, it would likely not suffice if the H<sub>2</sub> densities inward of 190 au were constrained and would likely become inconsistent when the height of the CS emission surface was taken into account. Therefore, the model from van Boekel et al. (2017) provides the most consistent profile in the outer disk. This is not surprising as this profile was found by fitting the radial profile of scattered light out to ~200 au, while the previous two surface density profiles were inferred from integrated flux values.

Although uncertainties in M*, H<sub>mid</sub>, H<sub>2</sub> ortho to para ratios and the location of the emission will propagate into the minimum mass, these are negligible compared to the lack in sensitivity of this approach to n<sub>H<sub>2</sub></sub> ≳ 10<sup>7</sup> cm<sup>-3</sup> due to the thermalization of the *J* = 7–6 transition. For commonly assumed power-law surface density profiles, 95% of the disk mass is within 80 au for TW Hya, far interior to where this technique is sensitive. Observations of higher density tracers, such as the *J* = 9–7 transition of CS would enable tighter constraints of Σ<sub>gas</sub> at a larger range of radii.

## 6. Summary

We have used spatially and spectrally resolved observations of the *J* = 7–6, 5–4, and 3–2 transitions of CS in TW Hya to constrain the physical conditions from where CS arises both in a spatially resolved manner and as a radial profile. Accounting for the rotation of the gas and azimuthally averaging the spectra allows us to apply the latter method out to 230 au.

Through both an LTE and a non-LTE approach to fitting the line ratios, we find an azimuthally symmetric physical structure. This approach demonstrated that the transitions were thermalized and thus provided a lower limit to the local H<sub>2</sub> density. The column density of CS was found to have two significant knees at 90 and 160 au; however, distinguishing between an abundance depletion of a true depletion in the total gas surface density cannot be done. We are able to place an upper limit on the nonthermal broadening component for all three lines in regions traced by the CS emission, *z* ≲ *H*, finding *M* ≲ 0.1 across the disk, consistent with previous determinations for this source (Hughes et al. 2011; Flaherty et al. 2018).

Extrapolating the H<sub>2</sub> volume density limits to a minimum gas surface density profile allowed us to place a limit on the minimum mass of the disk of 3 × 10<sup>-4</sup> M<sub>Sun</sub>, in line with constraints from molecular line emission including HD emission (Bergin et al. 2013; Trapman et al. 2017), in addition to constraining the gas surface density profile in the outer disk. Observations of higher *J* transitions will extend the sensitivity

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Continuum flux measurements for TW Hya at multiple wavelengths. The blue points show measurements from the data present in this paper while black are literature values.</figcaption></figure>

of this method to larger densities and thus allow for tighter constraints on the disk mass.

This paper serves as a template for future multiband observations and demonstrates the power of line excitation analyses in determining spatially resolved physical structures.

We thank the anonymous referee for comments that have improved the clarity of this paper. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00440.S and ADS/JAO.ALMA#2013.1.00387.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/ NRAO, and NAOJ. This work was supported by funding from NSF grants AST-1514670 and NASA NNX16AB48G. R.T. would like to thank Ryan Loomis for helpful discussions and Adam Ginsburg for his extensive help with pyradex.

## Appendix A
Flux Calibration

To check for any significant differences that may arise in the absolute flux scaling of the data, we compare the continuum flux values measured in each of our bands, 0.138, 0.603, and 1.430 Jy at 134, 242, and 353 GHZ, respectively, with previously published data. We plot these in Figure 8 using flux measurements from Qi et al. (2004, 2006), Hughes et al. (2009), Tsukagoshi et al. (2016), and Zhang et al. (2016).

No significant deviation is found relative to the literature values, suggesting that the flux calibration was performed adequately.

## Appendix B
Azimuthally Averaged Spectra

Figure 9 shows the azimuthally averaged spectra for the three transitions. Within the inner ≈25 au, the beam smear results in highly non-Gaussian line profiles; however, outside this radius a Gaussian profile is an apt description. We detect *J* = 3–2 emission out to 230 au, *J* = 5–4 out to 220 au, and *J* = 7–6 out to 210 au. These radii are consistent with the outer edge of <sup>12</sup>CO at 215 au (Huang et al. 2018). In Figure 10, we show the residuals from a Gaussian fit to the data. We calculate the noise at each radius taking into account the number of beams averaged over to produce the spectrum. While at smaller radii there appear to be large systematic deviations from a
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Azimuthally averaged spectra centered on $V_{LSR}$ = 2.842 km s⁻¹. The label on the right shows the radius. There is a 3 K offset between each radius.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Residuals from the azimuthally averaged spectra shown in Figure 9 and a Gaussian line profile. The shading behind each residual shows the noise for that radius taking into account the stacking of independent spectra.</figcaption></figure>

Gaussian profile due to the beam smearing, they are of comparable magnitude to the noise and should thus not significantly bias the results.

## Appendix C
### Posterior Distributions for Non-LTE Modeling

In this appendix, we discuss the MCMC sampling of the posterior distributions used in the non-LTE fitting described in Section 4. We take a representative selection at 170 au. Figures 11, 12, and 13 show the posterior distributions of the excitation conditions and noise models, respectively. At the top of each column of panels is the median value of the distribution with the uncertainties corresponding to the 16th and 84th percentiles of the distribution.

Figure 11 shows that both $T_{\text{kin}}$ and $N(\text{CS})$ are well constrained and only slightly correlated, while only limits can be placed on $n(\text{H}_2)$ and $\log_{10} M$. The steep fall-off of the $\log_{10} M$ posterior distribution demonstrates that a tight upper limit has been found.

For the noise properties, no correlation between the parameters are observed. $\sigma$, with units of kelvin, and demonstrate the significant increase in SNR achieved through the azimuthal averaging compared to the channel noise reported in Table 1. The correlation length, $\ell$, are in units of km s$^{-1}$, yielding $\ell/\Delta V$ values of between 3 and 5. This is consistent with the noise seen in Figure 9 where noise features appear to be correlated over 3–5 channels.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Posterior distributions of the excitation conditions from a representative fit at 170 au when fitting with $\log_{10} M$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12. Posterior distributions of the excitation conditions from a representative fit at 170 au when fitting with a linear *M*.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13. Posterior distributions of the noise models from a representative fit at 170 au.</figcaption></figure>

## ORCID iDs

Richard Teague https://orcid.org/0000-0003-1534-5186

Stéphane Guilloteau https://orcid.org/0000-0003-3773-1870

Edwin A. Bergin https://orcid.org/0000-0003-4179-6394

Mario Flock https://orcid.org/0000-0002-9298-3029

Uma Gorti https://orcid.org/0000-0002-3311-5918

## References

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL, 820, L40

Bai, X.-N. 2015, ApJ, 798, 84

Bai, X.-N. 2017, ApJ, 845, 75

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Balbus, S. A., & Hawley, J. F. 1998, RvMP, 70, 1

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Natur, 493, 644

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101

Boehler, Y., Weaver, E., Isella, A., et al. 2017, ApJ, 840, 60

Calvet, N., D'Alessio, P., Hartmann, L., et al. 2002, ApJ, 568, 1008

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204

Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496

Czekala, I., Mandel, K. S., Andrews, S. M., et al. 2017, ApJ, 840, 49

Denis-Alpizar, O., Stoecklin, T., Guilloteau, S., & Dutrey, A. 2018, MNRAS, 478, 1811

Dipierro, G., Ricci, L., Pérez, L., et al. 2018, MNRAS, 475, 5296

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130

Fang, T., & White, M. 2004, ApJL, 606, L9
Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJL, 776, L38

Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68

Flower, D. R., & Watt, G. D. 1984, MNRAS, 209, 25

Flower, D. R., & Watt, G. D. 1985, MNRAS, 213, 991

Foreman-Mackey, D., Agol, E., Angus, R., & Ambikasaran, S. 2017, AJ, 154, 220

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Forgan, D., Armitage, P. J., & Simon, J. B. 2012, MNRAS, 426, 2419

Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343

Gammie, C. F. 2001, ApJ, 553, 174

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Gorti, U., Hollenbach, D., Najita, J., & Pascucci, I. 2011, ApJ, 735, 90

Guilloteau, S., Dutrey, A., Wakelam, V., et al. 2012, A&A, 548, A70

Guilloteau, S., Piétu, V., Chapillon, E., et al. 2016, A&A, 586, L1

Hendler, N. P., Pinilla, P., Pascucci, I., et al. 2018, MNRAS, 475, L62

Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018, ApJ, 852, 122

Hughes, A. M., Wilner, D. J., Andrews, S. M., Qi, C., & Hogerheijde, M. R. 2011, ApJ, 727, 85

Hughes, A. M., Wilner, D. J., Cho, J., et al. 2009, ApJ, 704, 1204

Kastner, J. H., Zuckerman, B., Weintraub, D. A., & Forveille, T. 1997, Sci, 277, 67

Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869

Lesur, G. R. J., & Latter, H. 2016, MNRAS, 462, 4549

Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17

Loomis, R. A., Cleeves, L. I., Öberg, K. I., et al. 2018, ApJ, 859, 131

Lyra, W., & Klahr, H. 2011, A&A, 527, A138

Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3

Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87

Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017, ApJ, 842, 9

Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20

Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012, A&A, 541, A97

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Sci, 353, 1519

Pohl, A., Benisty, M., Pinilla, P., et al. 2017, ApJ, 850, 52

Qi, C., Ho, P. T. P., Wilner, D. J., et al. 2004, ApJL, 616, L11

Qi, C., Wilner, D. J., Calvet, N., et al. 2006, ApJL, 636, L157

Rohlfs, K., & Wilson, T. L. 1996, Tools of Radio Astronomy (New York: Springer-Verlag)

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013, ApJ, 775, 73

Simon, J. B., Hughes, A. M., Flaherty, K. M., Bai, X.-N., & Armitage, P. J. 2015, ApJ, 808, 180

Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49

Teague, R., Semenov, D., Gorti, U., et al. 2017, ApJ, 835, 228

Thi, W. F., Mathews, G., Ménard, F., et al. 2010, A&A, 518, L125

Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69

Tsukagoshi, T., Nomura, H., Muto, T., et al. 2016, ApJL, 829, L35

van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

Walsh, C., Daley, C., Facchini, S., & Juhász, A. 2017, A&A, 607, A114

Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204

Zhang, K., Bergin, E. A., Blake, G. A., et al. 2016, ApJL, 818, L16

Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, NatAs, 1, 0130