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

# A Surface Density Perturbation in the TW Hydrae Disk at 95 au Traced by Molecular Emission

R. Teague¹, D. Semenov¹, U. Gorti²,³, S. Guilloteau⁴,⁵, Th. Henning¹, T. Birnstiel¹, A. Dutrey⁴,⁵,
R. van Boekel¹, and E. Chapillon⁴,⁵,⁶

¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany; teague@mpia.de

²NASA Ames Research Center, Moffett Field, CA 94035, USA

³SETI Institute, Mountain View, CA 94043, USA

⁴University Bordeaux, LAB, UMR 5804, F-33270 Floirac, France

⁵CNRS, LAB, UMR 5804, F-33270 Floirac, France

⁶IRAM, 300 rue de la Piscine, Domaine Universitaire, F-38406 Saint Martin d'Hères, France

Received 2016 October 24; revised 2016 December 7; accepted 2016 December 21; published 2017 January 31

## Abstract

We present ALMA Cycle 2 observations at 0.′′5 resolution of TW Hya of CS J = 5 − 4 emission. The radial profile of the integrated line emission displays oscillatory features outward of 1.′′5 (≈90 au). A dip-like feature at 1.′′6 is coincident in location, depth, and width with features observed in dust scattered light at near-infrared wavelengths. Using a thermochemical model indicative of TW Hya, gas-grain chemical modeling, and non-LTE radiative transfer, we demonstrate that such a feature can be reproduced with a surface density depression, consistent with the modeling performed for scattered-light observations of TW Hya. We further demonstrate that a gap in the dust distribution and dust opacity only cannot reproduce the observed CS feature. The outer enhancement at 3.′′1 is identified as a region of intensified desorption due to enhanced penetration of the interstellar far-UV radiation at the exponential edge of the disk surface density, which intensifies the photochemical processing of gas and ices.

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

## 1. Introduction

Protoplanetary disks are sites of active planet formation, and it is in the dense midplanes of these disks where grain growth seeds planetesimal formation. Once a planet has grown to a certain mass, it will begin interacting with the disk and sculpt the density structure (Kley & Nelson 2012; Turner et al. 2014). The hallmarks of ongoing planet formation and planet–disk interactions are frequently observed, such as spirals and gaps in scattered light (Debes et al. 2013; Benisty et al. 2015; Wagner et al. 2015) or rings and dust traps seen in thermal millimeter-dust continuum (van der Marel et al. 2013; ALMA Partnership et al. 2015; Andrews et al. 2016).

While embedded planets are an attractive mechanism for producing these features, other physical mechanisms unrelated to planets have been shown to produce similar features, such as gravitational instabilities, dead zones of the magnetorotational instability (MRI), snowlines of volatile species or dust evolution (Birnstiel et al. 2015; Flock et al. 2015; Pohl et al. 2015; Zhang et al. 2015). To connect these disparate regions requires studies of molecular line emission, which traces the gas between the areas probed by thermal continuum or scattered light.

TW Hya, as the nearest protoplanetary disk to us at 59 ± 1 pc (Gaia Collaboration et al. 2016), is the most well studied disk across a huge wavelength range and is the ideal target for such a study. The disk displays a variety of substructures and radial morphologies in scattered light (Debes et al. 2013; Akiyama et al. 2015; Rapson et al. 2015; van Boekel et al. 2017), millimeter-dust continuum (Andrews et al. 2012, 2016; Menu et al. 2014; Tsukagoshi et al. 2016), and molecular line emission (Qi et al. 2013; Kastner et al. 2015; Bergin et al. 2016; Nomura et al. 2016; Schwarz et al. 2016), with many physical and chemical processes being invoked to explain these striking features.

In this paper we present ALMA Cycle 2 observations of CS J = 5 − 4 line emission, whose radial profile shows a shallow dip coincident with the ~90 au dust gap seen in the scattered light (Debes et al. 2013; Rapson et al. 2015; van Boekel et al. 2017). Section 2 describes the observations and the observed features, while in Section 3 we follow the modeling carried out by van Boekel et al. (2017) and argue that these features are the product of a depression in the surface density of the gas rather than a chemical effect. Our findings are discussed in more detail, including their power as complimentary studies for future high-resolution observations, in Section 5. A summary follows in Section 6.

## 2. Observations

### 2.1. Data Reduction

The observations targeted CO J = 2 − 1 at 230.538 GHz, CN N = 2 − 1 at 226.875 GHz, and CS J = 5 − 4 at 244.936 GHz and were originally presented in Teague et al. (2016). However, we provide a brief description of the data here and focus only on the CS J = 5 − 4 emission.

The ALMA Cycle 2 observations, project 2013.1.00387.S, were performed on 2015 May 13. Data were calibrated using the standard ALMA calibration script in the CASA software package.<fnref n="7" /> The absolute flux calibration was estimated using Ganymede. The derived flux for our amplitude and phase calibrator, J1037-2934, was 0.72 Jy at 228 GHz at the time of the observations, with a spectral index α = −0.54, while the ALMA flux archive indicated a flux of 0.72 ± 0.05 Jy between April 14 and 25. We hence estimate that the calibration uncertainty is about 7%.

<footnote n="7">⁷ http://casa.nrao.edu/</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Zeroth-moment map of CS J = 5 − 4 in TW Hya. The contours are in steps of 3σ, where σ = 5.35 mJy beam⁻¹ km s⁻¹, and the beam, with size 0′′ 42 × 0′′ 54, is shown in the bottom left corner. (b) Zeroth-moment map from panel (a) with an r² scaling applied to highlight the features in the outer disk. Such a scaling is used only to highlight the outer disk and is not physically motivated. A dip is clearly seen at r ≈ 1′′ 6, coincident with the dip in scattered-light emission (Debes et al. 2013; van Boekel et al. 2017).</figcaption></figure>

Figure 1 shows the continuum-subtracted zeroth moment of the CS J = (5 − 4) line emission in panel (a). No strong azimuthal structure is observed to be significant when considering the noise. The integrated flux of the line was 1.24 Jy beam⁻¹ km s⁻¹, calculated from an 8″ × 8″ box centered on the source. Panel (b) shows the same moment map, but with an r² scaling applied in order to highlight the outer regions of the disk. This scaling is applied to scattered-light data to account for the drop in stellar flux at a given radial position; however, for molecular line emission, such as the case here, this is not physically motivated and is applied purely to highlight the outer regions of the disk.

As no azimuthal structure is seen in the zeroth-moment
maps, nor suggested in previous observations, these data were
azimuthally averaged assuming an inclination of 7° (Qi
et al. 2004) and position angle 240° (Teague et al. 2016) to
yield a radial intensity profile, shown in Figure 2. The solid line
shows the zeroth-moment profile, while the dotted line shows
the r² scaling (panels (a) and (b) from Figure 1, respectively).
Gray shading around the lines demonstrates the 1σ noise in the
image. This does not take into account the estimated flux
calibration uncertainty of 7%, which would be a global
rescaling of the profile.

## 2.2. Radial Features

The peak of the emission is offset from the center at ≈1′′ 2.
Within this radius the continuum emission dominates the weaker
CS emission; thus, continuum subtraction may affect the assumed
inner radii of CS. Outside of this, the emission falls off and is
detected to ~3′′ 7 (Teague et al. 2016), comparable to that of CO
and scattered light (Andrews et al. 2012; Debes et al. 2013).

The integrated intensity emission profile shows oscillatory
features outward of 1′′5, as noted in Figure 2. Two features are
present, an inner and an outer feature, each consisting of a
“dip” and a “bump” depending on the sign of the second
derivative, as shown in Figure 2(b). As there is no baseline with

<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Radial profile of the CS J = 5 − 4 emission from Figure 1(a). Annotated are the potential dip and enhancement locations. The beam FWHM is shown in the top right corner. (b) Second derivative of the radial profile highlighting the features. The gray regions show 1σ standard deviation of each azimuthal bin.</figcaption></figure>

which to make a comparison, these features could be
considered dips in emission at 1"'6 and 2"'6, enhancements at
2"'2 and 3"'2, or a combination of both.

As the millimeter continuum has a sharp truncation at
~60 au (Hogerheijde et al. 2016), a comparison with the
substructures observed by Andrews et al. (2016) is impossible.
Scattered-light emission, however, extends out to ~220 au,
making for a better comparison (Akiyama et al. 2015; Rapson
et al. 2015; van Boekel et al. 2017). The most prominent
feature, “Dip I,” at 1"'6, or 95 au when assuming a distance of
d = 59 pc (Gaia Collaboration et al. 2016),<fnref n="8" /> matches the
location and apparent depth to the dip in scattered-light
emission previously observed with NICMOS on the Hubble
Space Telescope (Debes et al. 2013) and later confirmed with
ground-based observations. Such a similarity is suggestive of a

<footnote n="8">⁸ At the pre-Gaia distance of 54 pc (van Leeuwen 2007), this would be 86 au.</footnote>
common origin for these two features, a scenario that is explored further in Section 3.

An alternative scenario is that the features are enhancements in emission. Similar outer enhancements in molecular emission have been observed in several disks, for example, the secondary ring of DCO$^+$ in IM Lup (Öberg et al. 2015) and in CO isotopologues in AS 209 and the inner region of TW Hya (Huang et al. 2016; Schwarz et al. 2016). These works have all invoked a secondary desorption front, either photo or thermal, which releases volatiles from the ices, resulting in enhanced emission locally. However, these scenarios require a change in the local grain properties to produce a change in either UV penetration or thermal structure (Cleaves 2016). This prerequisite is satisfied for all three cases where the enhancement is found near the edge of the millimeter continuum, a proxy for the edge of the millimeter grains. For the two enhancements described here, “Bump I” and “Bump II,” the edge of the millimeter grains is too far inward to be the cause. However, “Bump II” is very close to the outer edge of the disk, and so may be due to enhanced desorption. This scenario is discussed in more detail in Section 5.

As much work has been done on secondary desorption fronts, in this paper we focus on the possibility of a feature induced by a surface density perturbation, as used in models of scattered-light observations.

### 3. Modeling Surface Density Perturbations

A commonly invoked mechanism to account for dips in scattered-light emission is a perturbation in the total surface density, a method used by both Debes et al. (2013) and van Boekel et al. (2017) to model the TW Hya scattered-light emission. In this section, we explore, using a model representative of TW Hya, the impact of such a surface density perturbation on the molecular emission of CS.

As molecular emission is a product of both excitation and abundance effects, it is important to perform self-consistent modeling taking into account the impact of the surface density perturbations on the physical structure and the resulting change in chemical abundances. To do this, we use an advanced thermochemical model in conjunction with chemical modeling and non-LTE radiative transfer to demonstrate the effect of a surface density perturbation on the radial profiles of the molecular emission and distinguish between scenarios.

#### 3.1. Physical Structure

For a baseline, we take the TW Hya model of Gorti et al. (2011) as the fiducial model. The unperturbed surface density is given by

$$ \Sigma_0 = 500 r^{-0.7} \exp\left(-\frac{r^{1.3}}{100}\right) \text{ g cm}^{-2}, \quad (1) $$

where r is the radial distance in au, yielding a total gas mass of $M_{\text{disk}} = 0.06 M_{\odot}$. The distribution was derived by fitting various gas spectra and dust continuum observed from the submillimeter to optical wavelengths. The model was made to reproduce integrated intensities rather than spatially resolved observations. Thus, while the model may not fully reproduce spatially resolved observations, it is a model with temperatures and densities expected in TW Hya. The model assumes a gas-to-dust ratio of 100, with an MRN-like grain size distribution, $n(a) \propto a^{-3.5}$, with maximum and minimum sizes of 1 mm and

0.09 µm, respectively. The assumed relative abundance of polycyclic aromatic hydrocarbons (PAH) per hydrogen nuclei is $10^{-9}$. The unperturbed model, Model A, is shown in Figure 3, with the gas density, gas temperature, and dust temperature shown in the panels, running left to right. The total mass of the unperturbed model is 0.05 $M_{\odot}$, consistent with the masses derived from HD measurements (Bergin et al. 2013).

The surface density is then perturbed with a Gaussian depression characterized by center $d_0$, width $\Delta d$, and depth $d$ such that the resulting surface density is given by

$$ \Sigma(r) = \Sigma_0(r) \times \left( 1 - d \cdot \exp\left[ -\frac{1}{2} \frac{(r - d_0)^2}{\Delta d^2} \right] \right). \quad (2) $$

We consider three perturbations to the surface density: Model B, a perturbation in the total surface density as in Debes et al. (2013), $d = 0.3$, $\Delta d = 20$ au, and $d_0 = 80$ au; Model C, the outer perturbation from van Boekel et al. (2017), $d = 0.55$, $\Delta d = 15$ au, and $d_0 = 80$ au; and Model D, a fully depleted gap, $d = 1$, $\Delta d = 9.45$ au, and $d_0 = 80$ au. The widths have been chosen such that each perturbed model has the same amount of mass removed relative to the unperturbed Model A: $2.7 \times 10^{-3} M_{\odot}$ (2.8 $M_{\text{Jup}$). This criterion for the width of the gap recovers the best-fit widths found with both Models B and C (Debes et al. 2013; van Boekel et al. 2017).

Given that the reduction of emission in scattered light could be due to a local depletion of small grains, we consider an additional set of models where the perturbation is applied only to the dust surface density. This results in the perturbed regions possessing an enhanced gas-to-dust ratio relative to the fiducial 100:1. With these models it is possible to explore the sensitivity of the molecular emission to the local available grain surface.

#### 3.2. Chemistry

With the prescribed surface densities covering between 3.9 and 200 au (TW Hya is observed to have material as close in as 1 au; however, this will not contribute to the low-energy molecular line emission we are concerned with here), the 1 + 1D disk physical structure was solved for self-consistently including heating and cooling processes following Gorti et al. (2011). Gas and dust temperatures are treated independently, allowing the gas temperatures in the upper, more strongly irradiated regions of the disk to deviate strongly from those of the dust. The stellar radiation was set up assuming a central stellar mass of $0.7 M_{\odot}$, a radius of $1.1 R_{\odot}$, and an effective temperature of 4200 K. This comprises a total far-UV (FUV) luminosity of $3 \times 10^{31} \text{ erg s}^{-1}$ and an X-ray spectrum that covers 0.1–10 keV and a total X-ray luminosity of $1.6 \times 10^{30} \text{ erg s}^{-1}$. The accretion rate is assumed to be $10^{-9} \text{ M}_{\odot} \text{ yr}^{-1}$, as constrained by the X-ray and UV observations (e.g., Brickhouse et al. 2012).

The resulting physical structures were used as the basis for the chemical modeling with the 1 + 1D code ALCHEMIC (Semenov et al. 2010), employing a full gas-grain network with deuterium fractionation, including up to triply deuterated species (Albertsson et al. 2014), and assuming uniformly sized amorphous olivine grains with the radius $a_d = 7$ µm. This grain size choice results in a grain surface per cell equivalent to that of the size ensemble used in the physical modeling.

In the chemical model, X-ray and UV photons and cosmic-ray particles (CRPs) are the sources of ionizing radiation. The standard CRP ionization rate of $1.3 \times 10^{-17} \text{ s}^{-1}$ is adopted
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Gas density, (b) temperature, and (c) dust temperature for the fiducial model, Model A.</figcaption></figure>

throughout the entire disk. The FUV penetration is computed in a 1 + 1D approximation as in Semenov & Wiebe (2011), using the visual extinction $A_V$ in the direction toward the central star for the stellar FUV component and the visual extinction $A_V^{IS}$ in the vertical direction for the interstellar (IS) FUV component. To compute visual extinction $A_V$, we used the following approximate scaling relation:

$$A_V (\text{mag}) = \frac{N(\text{H})(\text{cm}^{-2})}{1.59 \times 10^{21}} \frac{0.1}{a_d (\mu\text{m})} \frac{m_{d/g}}{10^{-2}}, \quad (3)$$

where $a_d$ is the grain radius in $\mu$m and $m_{d/g}$ is the dust-to-gas mass ratio. The unattenuated stellar FUV intensity at 100 au is $\chi_\star = 410\chi_0$, as expressed in the Draine (1978) FUV field units. The unattenuated interstellar FUV intensity $\chi_{IS} = 1\chi_0$. Then, the FUV flux at a given disk location could be roughly computed as

$$I_{\text{UV}} = \chi_{\star} e^{-A_{V}^{*}/1.086} + \chi_{\text{IS}} e^{-A_{V}^{\text{IS}}/1.086}. \quad (4)$$

To compute the X-ray ionization rate, we used the parameterization of Bai & Goodman (2009), as presented in Equation (37) from Armitage (2015) for the average X-ray photon energy of 3 keV.

A “low metals,” mainly atomic set of initial abundances was used (Lee et al. 1998; Semenov et al. 2010). Table 1 details the relative abundance of the initial species.

Chemistry was modeled over 1 Myr without taking dynamics and disk evolution into account, which is an appropriate assumption for such chemically simple, fast-evolving species as CO, CN, and CS (Semenov & Wiebe 2011).

### 3.3. Radiative Transfer

These chemical abundances were then used for non-local thermal equilibrium (non-LTE) 3D radiative transfer modeling using LIME (Brinch & Hogerheijde 2010). In the modeling we assumed a turbulent broadening component of $v_{turb} = 0.3c_s$ where $c_s$ is the local sound speed (Teague et al. 2016). The collisional rates for CS were taken from Lique et al. (2010), stored in the LAMDA database.<fnref n="9" />

A final step was to make a comparison with the observation. The simobserve task in CASA was used to simulate observations using the same array configuration in order to account for the uv sampling and spatial resolution.

<table><caption>**Table 1**
Initial Chemical Abundances</caption><thead><tr><th>Species</th><th>Abundance<sup>a</sup><br>(n(X)/n(H))</th><th>Species</th><th>Abundance<sup>a</sup><br>(n(X)/n(H))</th></tr></thead><tbody><tr><td>ortho-H<sub>2</sub></td><td>3.75 (-1)</td><td>S</td><td>9.14 (-8)</td></tr><tr><td>para-H<sub>2</sub></td><td>1.13 (-1)</td><td>Si</td><td>9.74 (-9)</td></tr><tr><td>He</td><td>9.75 (-2)</td><td>Fe</td><td>2.74 (-9)</td></tr><tr><td>O</td><td>1.80 (-4)</td><td>Na</td><td>2.25 (-9)</td></tr><tr><td>C</td><td>7.86 (-5)</td><td>Mg</td><td>1.09 (-8)</td></tr><tr><td>N</td><td>2.47 (-5)</td><td>Cl</td><td>1.00 (-9)</td></tr><tr><td>HD</td><td>1.55 (-5)</td><td>P</td><td>2.16 (-10)</td></tr></tbody></table>

**Note.**

<sup>a</sup> Read *a(b)* as *a* × 10<sup>-b</sup>.

## 4. Results

### 4.1. Physical Structure

A depression in the surface density will create a change in the vertical structure of the disk, altering the amount of incident radiation and consequently the temperature and density of the region, impacting the chemical evolution. Figure 4 compares the different $n(H_2)$ and $T_{gas}$ structures for the seven models around the perturbation center. As $T_{dust}$ remains largely unchanged, with average deviations from Model A of 2% for $z/r \le 0.75$, we do not plot those figures. A full set of physical structure plots can be found in the Appendix in Figures 10–16.

The shallow ($d = 0.3$) perturbation in Models B and Bd makes no significant change to the physical structures compared to the physical structure of the standard disk, Model A. A slight reduction of the surface density over the dip does not affect the disk temperature structure, while gas density variations over the dip are too weak to be visible in the top row of Figure 4.

With the slightly deeper perturbation in Model C, the scale height of the disk is reduced and results in a less vertically extended disk. With the depletion of smaller grains, the gas is unable to cool as efficiently in Model Cd relative to Model C, resulting in a slightly warmer region with an average temperature increase of ≈5 K at the center of the perturbation, 80 au, for $z/r \le 0.75$.

The effects of a fully depleted region in Models D and Dd are much more prominent, as shown in the final two columns of Figure 4. For Model D, the reduction in surface density lowers the

<footnote n="9">⁹ http://home-strw.leidenuni.nl/~moldata/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of the *n*(H₂) (top) and T<sub>gas</sub> (bottom) structures for the seven models over the gap location, 80 ± 30 au. The physical structure outside this radial region is comparable for all seven molecules. The dust temperature structure is indistinguishable between the models and so is not shown.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Column densities of CS from Model A in gray, compared to those from the perturbed models. Red lines are for models with the perturbation applied to the total surface density; blue lines are where the perturbation is just in the dust surface density. Models B, C, and D are in panels (a), (b), and (c), respectively.</figcaption></figure>

scale height of the disk, resulting in a highly irradiated outer surface. This enhances UV-driven processes, such as H₂ pumping and photoelectric and PAH heating, resulting in a local temperature increase. In the dust-depleted region in Model Dd, collisional cooling becomes less efficient, and gas gets locally hotter than in the unperturbed Model A. This results in a sharp spike in gas temperature, also increasing the gas scale height locally.

## 4.2. Chemical Structure

In this section we restrict our discussion to the column densities for CS. A more extensive exploration of the whole chemical inventory is left for future work. Figure 5 shows the column densities of CS for the seven models considered.

In the chemical network used, CS is efficiently formed through the recombination of HCS⁺, or HOCS⁺ at warmer temperatures, T > ∼ 50 K, with H₂. It is primarily depleted from the gas phase by freezing out onto grain surfaces. The predicted column density for our unperturbed Model A at r = 100 au is N(CS) ~ 10¹³ cm⁻² and is in good overall agreement with typically observed values in T Tauri disks (Guilloteau et al. 2012) and estimated for TW Hya from single-dish observations (Kastner et al. 2014).

Models B and Bd show little change from the fiducial model, an expected outcome given the lack of change in physical structure. This demonstrates that changes in volume density on the order of 30% are insufficient to impact the CS chemistry. In contrast, Model C shows a clear depletion in N(CS) relative to Model A. With the reduced dust density of the disk in Models
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Continuum-subtracted normalized radial profiles of CS J = 5 – 4 for the seven models (colored lines) compared to TW Hya (black line). The models have been processed using the *simobserve* task in CASA to match the conditions of the observations. A normalized profile is used as this modeling only aims to reproduce the radial morphology of the emission profile rather than the absolute intensity.</figcaption></figure>

C and Cd, the harsher stellar UV radiation dominates slightly deeper into the disk, truncating the vertical extent of the CS molecular layer, resulting in a narrower CS molecular layer.

The enhancement of Model Cd relative to Model C is due to the lack of grain surface onto freeze-out in the region with an enhanced gas-to-dust ratio. This reduced efficiency of depletion is more clearly seen in Model Dd, which shows a spike of N(CS) where there are no grains available for the freeze-out, in stark contrast to what is observed in TW Hya. Furthermore, a fully depleted region of gas results in no molecular gas, as shown in Model D with a large gap in N(CS). Such extreme changes in column density should be readily detectable in the emission profiles.

### 4.3. Molecular Emission

Figure 6 shows the CS J = 5 – 4 radial intensity profiles after continuum subtraction for all seven models. The left column shows the radial profile directly from LIME, and the right column shows the radial profile for the simulated observations. The solid lines represent models with constant gas-to-dust ratios, while dotted lines are for models with dust surface density perturbations, and the shaded regions show the azimuthal variance of the radial bin. For the LIME model, the noise is intrinsic Monte Carlo noise from the random grid in LIME.

The intensities in the inner disk are roughly a factor of 2 too low compared to the observations, although they are more comparable for the outer ($r \gtrsim 100$ au) disk. For this work, we focus on the deviation between models, rather than a direct fit of the TW Hya disk, and discuss the applicability of the model to TW Hya in Section 5.

As expected from the column densities, Models B and Bd exhibit no clear deviation from the unperturbed Model A, further demonstrating that a shallow dip cannot simultaneously explain the CS emission morphology and the NICMOS observations of Debes et al. (2013). Conversely, Model D and Dd produce clear features that would be easily recognizable in the observations had the surface density perturbation been so dramatic.

Models C and Cd show a clear deficit of emission relative to Model A due to the reduced column densities described in the

previous section. Both emission profiles here show a more comparable morphology to the observations than other models. This is more clearly seen in Figure 7, which compares the normalized radial profiles with the observations. Panel (a) shows the LIME models, and panel (b) shows the simulated observations. Model C clearly has a radial morphology that is comparable to the observations, while Model Cd lacks the depression at 1.6.

These simulated observations demonstrate that the surface density perturbation used in van Boekel et al. (2017) is able to produce a similar radial morphology in CS J = 5 – 4 emission to that observed in TW Hya. Furthermore, the more shallow perturbation used in Debes et al. (2013) was insufficient to make a noticeable change from the fiducial model, and a fully depleted gap would result in features that are far more prominent in the observations.

## 5. Discussion

In the previous section we have demonstrated that a surface density perturbation, which has been found to well model the scattered-light profile of TW Hya, also produces a similar feature in the CS emission that has been observed in TW Hya. In the following section we first discuss the applicability of the results to TW Hya before discussing other possible mechanisms for producing the oscillatory features we observe.

### 5.1. Applicability to TW Hya

Clearly these models do not fully recover the observation profiles of TW Hya: the CS J = 5 – 4 line has a peak flux of 34.5 mJy beam⁻¹ km s⁻¹ at a radial offset 0.58, while the fiducial model, Model A, has a peak flux of 26.7 mJy beam⁻¹ km s⁻¹ at 1.01. While an iterative modeling effort was not the goal of this work, it is important to understand what may cause the discrepancy given the “TW Hya”-like fiducial model and whether this would impact the conclusions.

To verify that the physical and chemical models were producing comparable flux densities to TW Hya, radiative transfer modeling was performed for ¹²CO, ¹³CO, and C¹⁸O in Model A. Abundances for isotopologues were scaled assuming n(¹²CO)/n(¹³CO) = 84, consistent with the interstellar
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparing the second derivative of the CC $J = 5 - 4$ normalized radial intensity profiles for the seven models (colored lines) with the observations (black line). This demonstrates that Model C produces the most consistent radial morphology with the observations.</figcaption></figure>

medium (Wilson & Rood 1994), and $n(^{13}\text{CO})/n(^{18}\text{O}) = 10$, which was found for TW Hya (Schwarz et al. 2016). Flux densities were calculated for the CO $J = 2 - 1, ^{13}\text{CO} J = 3 - 2, 6-5$, and $^{18}\text{O} J = 3 - 2, 6-5$ lines, which all agreed within 10% of the observed values (K. Schwarz 2017, private communication). This gives confidence that the model used is yielding realistic abundances of the targeted molecules, despite not accurately modeling TW Hya.

We note that while the modeled lines of CO isotopologues yield a comparable flux density, all fail to match the observed radial profile in the same manner: emission in the outer disk ($r \gtrsim 90$ au) is overproduced, while in the inner disk it is underproduced. This is likely because the fiducial model was made to fit integrated intensities, rather than spatially resolved observations, so that the exact radial profiles in temperature and density are not matched. Nonetheless, the modeling in Section 3 demonstrates that by perturbing an otherwise smooth surface density profile, features in the CS emission profile can be produced that are comparable in morphology to those that are observed in TW Hya.

## 5.2. An Outer Enhancement in CS

The observations also display a feature farther out in the disk between 2′′5 and 3′′5, “Dip II” and “Bump II.” This feature is reproduced in all models to varying extent, as shown in Figure 7. Given that all models reproduce this feature, it is likely intrinsic to the setup of the models and is attributed to the exponential edge of the surface density profile used (Equation (1)). With the surface density dropping precipitously outside 100 au, UV radiation is able to penetrate deeper into the disk, increasing the desorption rate of volatiles.

Figure 8 shows the local UV flux for Model A with the CS molecular layer, defined as cells that contribute 1% or more to the CS column density at that radius, overlaid. Although the molecular layer traces a lower UV flux in the outer disk, the CS number density weighted average UV flux increases radially outward from 140 au for all models. This relative enhancement

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Local UV flux within Model A in Habing units. The hatched region shows the CS molecular layer, defined as cells that contribute 1% or more to the total column density. The white dotted contours show $n(\text{H})$ in cm$^{-3}$.</figcaption></figure>

in UV flux increases the intensity of photoprocessing, which impacts the CS abundance.

While photodesorption of CS could directly enhance the CS abundance, other pathways are possible. Enhanced UV intensities could lead to a higher rate of dissociation of CO, releasing more C with which to form light hydrocarbons, in addition to unlocking elemental S from other more complex species. Furthermore, higher ionization rates from the increased FUV intensity could ionize S to form $S^+$, resulting in much more efficient formation of CS.

While for the chemical models we can isolate the enhanced UV intensity as the cause, potentially resulting in a photo-desorption front, a similar mechanism to that proposed by Öberg et al. (2015), we cannot make such a claim for TW Hya. An outer enhancement could equally be due to increased thermal desorption front (e.g., Cleeves 2016). To further our understanding of this feature requires a much better characterization of the outer disk. This can be achieved with much more sensitive observations of more complex species than CO isotopologues that only extend out to ~90 au (Schwarz et al. 2016).
### 5.3. Embedded Planet

An attractive and commonly invoked mechanisms for explaining surface density perturbations is a planet-carved gap (Lin & Papaloizou 1986). Crida et al. (2006) show that the depth and width of a depression in the surface density can be used to constrain the mass of embedded planets. Debes et al. (2013) used this relationship with their perturbed surface density to estimate an embedded planetary mass of between 6 and 28 $M_⊕$.

For the partially filled case, as suggested by our modeling, Duffell (2015) derives an analytic expression for the planet mass,

$$q^2 = \frac{3\pi\alpha d}{(1-d)f_0 M^5}, \quad (5)$$

where $q = m_{\text{planet}}/M_*$, $\alpha$ is the viscosity parameter (Shakura & Sunyaev 1973), $d$ is the depth as in Equation (2), $f_0 = 0.45$ (a dimensionless parameter), and $M = r/H_p$ is the Mach number. Taking the depth of the perturbation from Model C, $d = 0.45$, consistent with the depth found by van Boekel et al. (2017), $M_* = 0.7 M_⊕$, and $\alpha = 10^{-4}-10^{-3}$, values appropriate for TW Hya (Teague et al. 2016), we find, using Equation (5), a potential planet mass of 12–38 $M_⊕$.

A better estimation could be made with a more precisely determined depth, $d$, requiring the dip feature to be observed in multiple molecular species equally sensitive to changes in local density and constraints on the local value of $\alpha$. In addition, high-resolution observations of HCN and its isotopologues have been proposed as chemical tracers of such an embedded planet (Cleeves et al. 2015a), allowing for a direct indication of a planet-opened gap.

### 5.4. Disk Instabilities

Pure hydrodynamical or magnetohydrodynamical instabilities can also create significant perturbations in the gas density structure (see for a review Turner et al. 2014). These instabilities can be broadly split into two types: ones producing azimuthally symmetric features, and those producing distinct azimuthal structures.

Molecular emission from TW Hya, in addition to the millimeter-continuum and scattered-light observations, shows no strong azimuthal variance (Debes et al. 2013; Qi et al. 2013; Akiyama et al. 2015; Kastner et al. 2015; Rapson et al. 2015; Andrews et al. 2016; Nomura et al. 2016; Schwarz et al. 2016; Tsukagoshi et al. 2016; van Boekel et al. 2017). Thus, instabilities that routinely produce vortices, spiral arms, or azimuthally asymmetric gaps are likely not present at ~90 au in TW Hya. These include the baroclinic instability (Klahr & Bodenheimer 2003; Cossins et al. 2010; Douglas et al. 2013; Takahashi et al. 2016), vertical shear instability (Nelson et al. 2013; Stoll & Kley 2014; Richard et al. 2016), zombie vortex instability (Marcus et al. 2015), and Rossby wave instability (Lovelace et al. 1999; Varnière & Tagger 2006; Lyra et al. 2015).

A more promising candidate is the MRI (Balbus & Hawley 1990). Flock et al. (2015) simulated a generic T Tauri disk assuming a gas-to-dust ratio of 100 that developed a surface density depression outside the dead zone with a depth of $d \approx 0.5$. The change in density structure and Elsässer number over this boundary forms axisymmetric zonal flows, carving the gap. A measurement of the radial extent of the dead zone (or

<figure><img src="image_10.png" /><figcaption>Figure 9. Radial intensity profiles for CS J = 5 – 4 from Model D assuming different beam sizes. All models have the same integrated intensity. An arbitrary offset has been included to allow for easier comparison between profiles.</figcaption></figure>

confirmation of its presence) would provide a great test for this instability; however, previous analyses of the ionization and turbulent velocity structures have yet to detect the edge (Cleeves et al. 2015b; Teague et al. 2016). Furthermore, Flock et al. (2015) note that no structure is seen in simulated scattered-light images, contrary to the case for TW Hya.

### 5.5. Grain Evolution

Motivated by the observed rings in scattered light, Birnstiel et al. (2015) demonstrated that grain evolution naturally leads to a region with a deficit of micron-sized grains without the need to invoke perturbed surface densities. This arises because of the lack of replenishment of micron-sized grains due to the slow growth and drift timescales at this region. Such a deficit of micron-sized grains manifests as a dip in scattered-light emission, but not thermal millimeter continuum, because of the lack of grain surface to reflect the light.

Such a feature would have implications for molecular abundances; however, detailed chemical modeling of such a mechanism, requiring the transport of frozen-out volatiles, is beyond the scope of this work. Modeling of recent observations of $C_2H$ and $C_3H_2$ in TW Hya (Bergin et al. 2016) has shown that treating grain evolution can lead to significant changes in the local abundances of volatile species and is therefore an important avenue of future exploration.

### 5.6. Molecular Tracers

Despite the large archive of molecular line observations in TW Hya, the radial feature has thus far only been observed in CS. $N_2H^+$ observations do not betray structure in the outer disk (Qi et al. 2013). However, the stark inner hole has been used as a proxy of the CO snowline at $\approx 30$ au, which was attributed as the cause of the inner depression of $^{13}$CO and $C^{18}$O observed by Schwarz et al. (2016) and also seen by Nomura et al. (2016). The CO isotopologue emission also showed an enhancement at $\approx 60$ au, consistent with a secondary thermal desorption front (Cleeves 2016). This suggests that an enhancement in CS at $\approx 120$ au due to a secondary desorption front is an unlikely scenario because of the large distance from any clear transition in grain properties.

Kastner et al. (2015) observed $C_2H$ emission in a ring between 45 and 120 au with the Sub-Millimeter Array, arguing that the
observed line intensities and radial location suggest a vertical segregation of grain sizes and considerable dust grain processing. Bergin et al. (2016) further argued with observations of C₂H and C₃H₂ that hydrocarbon rings require a local enhancement of the C/O ratio, which could be achieved by the sedimentation and radial drift of grains rich in volatile ices. While these observations are not suggestive of a surface density perturbation, they do indicate that chemical abundances and the grain evolution are tightly coupled and must be considered together.

The modeling performed in Section 3 suggested that molecular lines such as N₂H⁺, H₂CO, HCN, and HNC are all sensitive to small changes in the volume density, as invoked here. A full treatment of the radiative transfer for these molecules is left for future work; however, it is clear that a larger repertoire of molecular tracers will provide unparalleled constraints on the physical structure of the outer disk.

### 5.7. Perturbation Depth

Better discrimination between physical mechanisms can be made with a more thorough characterization of the physical conditions over the gap region. As discussed above, this can be achieved through observations of molecular line emission. However, tracing features in molecular line emission rather than in continuum or scattered light will always be limited by the spatial resolution achievable. To understand how effective large beams are at masking such radial features, simulated observations of Model D were produced using a range of beam sizes, $\theta_{\text{beam}} = 0''2, 0''3, 0''5, \text{ and } 1''$, as shown in Figure 9.

This demonstrates that a feature is only detectable by visual inspection for beam sizes of comparable size or smaller than the feature size. Thus, with a gap FWHM of ≈0''4 at 59 pc, only observations with spatial resolutions of $\theta_{\text{beam}} \lesssim 0''5$ would potentially detect the feature, while a partially depleted gap, such as Model C, would need to be observed at an even higher resolution in order to detect the feature. Thus, not only is it essential for a range of molecular species to be observed, but also high-resolution observations are needed to resolve the features that may be hidden with large beam sizes.

## 6. Conclusion

We have reported observations of CS $J = 5-4$ emission from the TW Hya disk, which exhibits oscillatory features outward of $1''5$ (≈90 au). A dip-like feature is observed at $1''6$, comparable in morphology and position to a feature observed in the NICMOS total intensity scattered-light observations (Debes et al. 2013) and SPHERE polarized intensities (van Boekel et al. 2017).

Through self-consistent thermochemical modeling and extensive chemical and radiative transfer modeling, we have demonstrated that a depression in surface density results in comparable morphology of the radial profile observed in CS at this location. A surface density perturbation with a depth of 55% relative to a fiducial model, as used by van Boekel et al. (2017) to model SPHERE observations of TW Hya, produces the best fit to the morphology of the radial profile. A less severe perturbation of 30%, used by Debes et al. (2013) to model NICMOS observations, was insufficient to strongly alter the

physical and chemical structure of the disk, while a fully depleted region produced features that were far too deep. Models with the perturbation applied only to the dust surface density, therefore enhancing the local gas-to-dust ratio, were shown to impede the freeze-out of CS, resulting in an enhancement of CS emission, in clear contrast to the observations.

In addition, these models, using an exponential edge to the surface density profile, naturally reproduce the outer enhancement at $3''$. Due to the lower densities in the outer disk, interstellar UV can penetrate deeper into the disk. This can bolster the CS abundance through either enhanced desorption of CS or the unlocking of more C and S from more complex species to form CS.

These models demonstrate that both physical and chemical effects can work in tandem to produce observable features in molecular emission profiles. Only by better characterizing the physical structure of a protoplanetary disk can we begin to accurately model the physical and chemical processes at play. This can be achieved by exploiting the range of molecular tracers observable with ALMA and combing these data with auxiliary data including scattered-light and continuum images.

We are thankful to the anonymous referee, whose comments greatly improved the quality of this manuscript. R.T. is a member of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg, Germany. T.B. acknowledges support from the DFG through SPP 1833 Building a Habitable Earth (KL 1469/13-1). D.S. acknowledges support from the Heidelberg Institute of Theoretical Studies for the project “Chemical Kinetics Models and Visualization Tools: Bridging Biology and Astronomy.” This research made use of NASA’s Astrophysics Data System. This paper makes use of the following ALMA data: 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 the National Programs PCMI and PNPS from INSU-CNRS. U.G. acknowledges support from NASA ADAP award NNX14AR91G. Resources supporting this work were provided by the NASA High-End Computing program through the NASA Advanced Supercomputing Division at Ames Research Center.

## Appendix
### Additional Figures

In Figures 10–16 we display the physical structures of the models used in this work: Models A, B, Bd, C, Cd, D, and Dd, respectively. In all figures the left panel shows the gas density, the central panel shows the gas temperature, and the right panel shows the dust temperature. These physical models were used as the basis for chemical modeling and radiative transfer models.
<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Gas density, (b) temperature, and (c) dust temperature for the fiducial model, Model A.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. (a) Gas density, (b) temperature, and (c) dust temperature for Model B.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. (a) Gas density, (b) temperature, and (c) dust temperature for Model Bd.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** (a) Gas density, (b) temperature, and (c) dust temperature for Model C.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** (a) Gas density, (b) temperature, and (c) dust temperature for Model C.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** (a) Gas density, (b) temperature, and (c) dust temperature for Model D.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16. (a) Gas density, (b) temperature, and (c) dust temperature for Model Dd.</figcaption></figure>

References

Akiyama, E., Muto, T., Kusakabe, N., et al. 2015, ApJ, 802, L17

Albertsson, T., Semenov, D., & Henning, T. 2014, ApJ, 784, 39

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

Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162

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

Armitage, P. J. 2015, arXiv:1509.06382

Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737

Balbus, S. A., & Hawley, J. F. 1990, BAAS, 22, 1209

Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6

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

Birnstiel, T., Andrews, S. M., Pinilla, P., & Kama, M. 2015, ApJ, 813, L14

Brickhouse, N. S., Cranmer, S. R., Dupree, A. K., et al. 2012, ApJ, 760, L21

Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25

Cleeves, L. I. 2016, ApJ, 816, L21

Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015a, ApJ, 807, 2

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

Cossins, P., Lodato, G., & Testi, L. 2010, MNRAS, 407, 181

Crida, A., Morbidelli, A., & Massett, F. 2006, Icar, 181, 587

Debes, J. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45

Douglas, T. A., Caselli, P., Ilee, J. D., et al. 2013, MNRAS, 433, 2064

Draine, B. T. 1978, ApJS, 36, 595

Duffell, P. C. 2015, ApJL, 807, L11

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

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2

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

Hogerheijde, M. R., Bekkers, D., Pinilla, P., et al. 2016, A&A, 586, A99

Huang, J., Öberg, K. I., & Andrews, S. M. 2016, ApJL, 823, L18

Kastner, J. H., Hily-Blant, P., Rodriguez, D. R., Punzi, K., & Forveille, T. 2014, ApJ, 793, 55

Kastner, J. H., Qi, C., Gorti, U., et al. 2015, ApJ, 806, 75

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

Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211

Lee, H.-H., Roueff, E., Pineau des Forets, G., et al. 1998, A&A, 334, 1047

Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846

Lique, F., Spielfiedel, A., Feautrier, N., et al. 2010, JChPh, 132, 024303

Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805

Lyra, W., Turner, N. J., & McNally, C. P. 2015, A&A, 574, A10

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

Menu, J., van Boekel, R., Henning, T., et al. 2014, A&A, 564, A93

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

Nomura, H., Tsukagoshi, T., Kawabe, R., et al. 2016, ApJL, 819, L7

Öberg, K. I., Furuya, K., Loomis, R., et al. 2015, ApJ, 810, 112

Pohl, A., Pinilla, P., Benisty, M., et al. 2015, MNRAS, 453, 1768

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

Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Sci, 341, 630

Rapson, V. A., Kastner, J. H., Millar-Blanchaer, M. A., & Dong, R. 2015, ApJL, 815, L26

Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571

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

Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42

Semenov, D., & Wiebe, D. 2011, ApJS, 196, 25

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

Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77

Takahashi, S. Z., Tsukamoto, Y., & Inutsuka, S. 2016, MNRAS, 458, 3597

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

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

Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. of Arizona Press), 411

van Boekel, R., Henning, T., Menu, J. & SPHERE Collaboration 2017, ApJ, in press (arXiv:1610.08939)

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Sci, 340, 1199

van Leeuwen, F. 2007, A&A, 474, 653

Varnière, P., & Tagger, M. 2006, A&A, 446, L13

Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJL, 813, L2

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJL, 806, L7