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

HAL
open science

# Gas and Dust Shadows in the TW Hydrae Disk

Richard Teague, Jaehan Bae, Myriam Benisty, Sean M. Andrews, Stefano Facchini, Jane Huang, David Wilner

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

## To cite this version:

Richard Teague, Jaehan Bae, Myriam Benisty, Sean M. Andrews, Stefano Facchini, et al. Gas and Dust Shadows in the TW Hydrae Disk. The Astrophysical Journal, 2022, 930, 10.3847/1538-4357/ac67a3. insu-03705349

HAL Id: insu-03705349

https://insu.hal.science/insu-03705349

Submitted on 27 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Gas and Dust Shadows in the TW Hydrae Disk

Richard Teague¹, Jaehan Bae², Myriam Benisty³,⁴, Sean M. Andrews¹, Stefano Facchini⁵, Jane Huang⁶,<fnref n="7" /> and David Wilner¹

¹Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA; richard.d.teague@cfa.harvard.edu

²Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

³Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁴Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

⁵Università degli Studi di Milano, Via Celoria 16, I-20133 Milano, Italy

⁶Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA

Received 2022 February 16; revised 2022 April 13; accepted 2022 April 13; published 2022 May 13

## Abstract

We present new observations of CO J = 2 − 1 emission from the protoplanetary disk around TW Hya. Emission is detected out to 240 au (4") and found to exhibit azimuthal variations up to 20% beyond 180 au (3"), with the west side of the disk brighter than the east. This asymmetry is interpreted as tracing the shadow previously seen in scattered light. A re-analysis of the multi-epoch observations of the dust shadow in scattered light from Debes et al. suggests that an oscillatory motion would provide a better model of the temporal evolution of the dust shadow rather than orbital motion. Both models predict an angular offset between the dust shadow and the gas shadow of up to ~100°. We attribute this offset to the finite rate at which dust grains and gas molecules can exchange heat, dominated by the collisional rate between gas molecules and dust grains, $t_{\text{coll}}$. The angular offsets derived are equivalent to collisional timescales that range from the near-instantaneous up to $t_{\text{coll}} \sim 10$ yr, depending on whether a straight or curved dust shadow, as suggested by Hubble Space Telescope observations reported by Debes et al., is adopted. The inferred range of $t_{\text{coll}}$ are consistent with those predictions based on representative gas densities, temperatures, gas-to-dust ratios and grain sizes. These results represent the first time empirical constraints can be placed on $t_{\text{coll}}$.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Interferometry (808); Planet formation (1241)*

## 1. Introduction

TW Hya is one of the most well studied protoplanetary disks, owing primarily to its unique combination of proximity to Earth (d = 60.1 pc; Gaia Collaboration et al. 2018) and near-face-on orientation (i ≈ 5°; Teague et al. 2019b). As such, it is often used as a young analog of the solar system, guiding our understanding of how planetary systems like our own may have formed.

For TW Hya, the radial distribution of dust grains has been exquisitely mapped, both at (sub-)millimeter wavelengths, tracing the thermal emission of the millimeter-sized grains at the disk midplane (Andrews et al. 2016; Huang et al. 2018; Macías et al. 2021), and in the optical and NIR, tracing the scattering from sub-micrometer grains high in the atmosphere of the disk (Krist et al. 2000; Weinberger et al. 2002; Apai et al. 2004; Roberge et al. 2005; Debes et al. 2013, 2016; Akiyama et al. 2015; Rapson et al. 2015; van Boekel et al. 2017). These have revealed a striking level of substructure—concentric gaps and rings on astronomical unit scales punctuate an otherwise smooth background distribution across the entire radial extent of the disk—both along the disk midplane as well as in the atmosphere.

In concert, sub-millimeter interferometers have facilitated a comprehensive study of the gas component of the disk. Physical properties such as the gas temperature structure (Calahan et al. 2021), column density (Gorti et al. 2011; Bergin et al. 2013; Schwarz et al. 2016; Zhang et al. 2017), and radial extent (Huang et al. 2018) have all been probed through a combination of empirical analyses and complex thermochemical modeling. Similarly, spectral surveys have revealed a stunning level of chemical complexity, yielding unique insights into the ionization structure of the disk (Cleeves et al. 2015), the location of potentially planet-forming snowlines (Qi et al. 2013), and the astrochemical heritage of solar system bodies (Walsh et al. 2016; Loomis et al. 2018; Canta et al. 2021).

With the data available for TW Hya continuing to grow, a new frontier has opened up: temporal variability. Using archival Hubble Space Telescope (HST) data, Debes et al. (2017) demonstrated that a previously detected shadow in the outer disk, $r \gtrsim 50$ au (Roberge et al. 2005), appears to move across the face of the disk at a rate of 22°.7 yr⁻¹ in an anticlockwise direction (increasing PA). Curiously, this rotation appears to be counter to the disk rotation direction inferred through the winding direction of spirals, which are assumed to be trailing, reported in Teague et al. (2019b). Recent simulations presented by Nealon et al. (2019, 2020) demonstrated that a precessing inner disk could also be a plausible scenario to explain the counter-rotation.

In this paper we present new observations of CO J = 2 − 1 emission from TW Hya at a 0.′′15 angular resolution which we describe in Section 2. In Section 3, we show that there are azimuthal variations, consistent with the previously detected dust shadows. Motivated by the location of the gas shadows, a

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

<footnote n="7">⁷NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
re-analysis of the reported movement of the HST shadows in Section 4 suggests an oscillatory motion is also a plausible model for the temporal evolution of the dust shadow. With this model in hand, Section 5 uses the relative offset between the gas and dust shadows to place unique constraints on the thermal coupling between gas and dust in the outer disk. These new findings are discussed in Section 6 and summarized in Section 7.

## 2. Observations

Observations were acquired as part of project 2018.A.00021. S (PI: Teague), designed for high spatial and spectral resolution observation of ¹²CO J = 2 − 1, CN N = 2 − 1, and CS J = 5 − 4. The correlator was set up to cover the ¹²CO and CS line with a channel spacing of 15.3 kHz (resulting in a 30.5 kHz resolution, equivalent to 40 m s⁻¹ and 230.538 GHz), while CN was observed at a coarser 61 kHz channel spacing (yielding a 122 kHz resolution, equivalent to 160 m s⁻¹) in order to have bandwidth sufficiently wide to fully cover the J = 5/2 − 3/2 and J = 3/2 − 1/2 fine-structure groups. A single spectral window in FDM mode centered at 241.5 GHz was used for continuum emission in order to self-calibrate the data. A serendipitous detection of CH₂CN in the continuum window was presented in Canta et al. (2021). In this paper we focus on the CO emission, leaving the analysis of CS and CN for future work.

### 2.1. Calibration

Short-spacing data (baselines of 15–500 m) were taken on 2019 April 4 with two executions including 47.1 min of on-source time each. During these executions the quasars J1037-2934 and J1147-3812 were observed for calibration purposes, the former acting as bandpass and flux calibrator and the latter as phase calibrator. Long-baseline data (baselines ranging between 15 m and 2.62 km) were taken on 2019 September 29. Only four out of the requested 12 executions were observed due to scheduling constraints. Each execution included 50.1 min of on-source time. During these observations J1037-2934 was observed for bandpass calibration and flux calibration and J1126-3828 for phase calibration. The total on-source time was 4.9 hr.

Initial calibration was performed using the standard pipeline procedure in CASA v5.8.0. The data were then self-calibrated following the procedure used in the DSHARP program (Andrews et al. 2018). In brief, all spectral windows were used, masking out any lines in each spectral window. These line-free observations were used to derive phase solutions which were then applied to the entire data set. Prior to combining the different executions, each execution was aligned to a common phase center and the continuum fluxes were compared. All executions yielded fluxes that were within 2% of one another, except for the final long-baseline execution which deviated by about 10%. This execution was rescaled using the gaincal task such that the total flux matched that in the first long-baseline execution. The continuum was subtracted using the uvcontsub task.

### 2.2. Imaging

Briggs weighting and a channel spacing of 40 m s⁻¹ (the spectral resolution of the data) were chosen for the imaging. A variety of robust parameters were considered, and a value of 0.5 was found to give a good trade-off between angular

resolution and high-quality imaging. This resulted in a synthesized beam of 0.′′19 × 0.′′17 (89°). A Keplerian mask was generated,<fnref n="8" /> making sure that all disk emission was contained within the mask. Following Teague et al. (2019b), a stellar mass of M$_* = 0.81 M⊙ and a viewing geometry described by i = 5°.8 and PA = 151° were adopted for this mask. After the imaging, a correction was applied to the image to account for the non-Gaussian synthesized beams due to the combination of several different array configurations (see the discussion in Czekala et al. 2021), as proposed by Jorsater & van Moorsel (1995). The resulting rms in a line-free channel was measured to be 0.6 mJy beam⁻¹. Integrating over the Keplerian mask we recover an integrated intensities of 18.31 ± 0.03 Jy km s⁻¹. This uncertainty does not include a systematic uncertainty of ~10% associated with the flux calibration of the data.

The package bettermoments<fnref n="9" /> was used to calculate the zeroth moment map of the data. This was calculated using the Keplerian CLEAN mask without applying any σ-clipping. The resulting moment map is shown in Figure 1(b). The analysis of additional moment maps is the focus of a companion paper (R. Teague et al., 2022, in preparation).

## 3. Gas Shadows

No clear azimuthal structure is seen in the zeroth-moment map. Thus, to reveal underlying substructure, an azimuthally averaged radial profile was calculated using GoFish (Teague 2019a), adopting the same geometrical properties as used for the CLEAN mask to deproject the data into concentric annuli. To determine the disk center, a Keplerian rotation pattern was fit to the rotation map using eddy (Teague 2019b). This resulted in small offsets relative to the image center of less than a pixel in size: ΔR.A. = 20.0 ± 0.1 mas and Δdecl. = 2.0 ± 0.1 mas. This offset was necessary to include as the centering of the image during the data reduction was performed by a Gaussian fit to the continuum emission. The fitting of a Keplerian rotation pattern allows for a more precise measure of the disk center as both a far larger disk area is used for the fit and the strong azimuthal dependence of the velocity pattern aids in breaking rotational degeneracies of the model. The data were binned into annuli with a width of 1/4 of the beam FWHM (≈50 mas) to limit the underlying gradient in integrated intensity from biasing the measurement, with the caveat that measurements will be spatially correlated on scales comparable to the beam FWHM. The resulting azimuthally averaged radial profile is shown in Figure 1(a), with a morphology consistent with that found for the *J = 3 − 2* line presented by Huang et al. (2018).

This radial profile was used to generate an azimuthally symmetric background model which was subtracted from the integrated intensity map shown in Figure 1(b). The residuals are shown in Figure 1(c). A clear east/west asymmetry is seen, with the dashed lines marking a rough boundary between the two regimes (the lines trace a PA = 39°). On top of this large-scale asymmetry, there is much substructure seen within 3′′, likely associated with the perturbations reported in the gas temperature and dynamics by Teague et al. (2019b). To confirm that this asymmetry is unrelated to a mis-specified center, the same residuals were calculated for a range of offsets

<footnote n="8">⁸ github.com/richteague/keplerian_mask</footnote>

<footnote n="9">⁹ github.com/richteague/bettermoments</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 1. Summary of the CO integrated intensity. (a) Azimuthally averaged radial profile of the zeroth moment map for the $J=2-1$ emission in black, and the $J=3-2$ emission from Huang et al. (2018) in dashed gray. The shaded region about each line shows the standard deviation of each annulus and the error bar in the lower left showing the beam FWHM for each data set. (b) $J=2-1$ moment map with the beam size shown by the ellipse in the lower left corner. The solid contours show a signal-to-noise ratio of $5 \times 2^n$ with $n = \{0, 1, \dots, 7\}$. (c) Residuals after subtracting the azimuthally averaged radial profile from the data. The color scaling has been adjusted such that residuals smaller than $1\sigma$ are white. Solid blue and red lines show negative and positive $2\sigma$ and $5\sigma$ contours. The dashed line is a guide to emphasize the global east/west asymmetry.</figcaption></figure>

in both the x- and y-direction of $\pm 0.5$ pix (greater than a factor of 10 larger than the statistical uncertainties on the derived source center). In all cases the east/west asymmetry is present.

The global asymmetry is reminiscent of the optical and NIR shadows observed with the HST and reported by Debes et al. (2017). To make a comparison, we followed the procedure in Debes et al. (2017) to determine the azimuthal location of the CO shadow. In short, we fit the functional form to each radial bin of the integrated intensity map,

$$I(r, \phi) = \delta I(r) \cdot \cos(\phi - \phi_{s,gas}(r) - \pi) + \langle I(r) \rangle, \quad (1)$$

where $\langle I(r) \rangle$ is the azimuthally averaged integrated intensity, $\delta I(r)$ is the radially dependent amplitude of the azimuthal variation in integrated intensity and $\phi_{s,gas}(r)$ is the radially dependent position angle of the shadow (note that unlike Debes et al. 2017 we have included a $\pi$ phase offset such that $\phi_{s,gas}$ describes the azimuthal minimum, not the peak). This form makes the assumption that the shadowing is smooth and affects the whole azimuth, rather than a narrow region with sharp, well defined edges (e.g., Casassus et al. 2019). Similarly, it should be noted that $\phi_{s,gas}$ does not measure the azimuthal minimum of each annulus, but rather the azimuthal location of the sinusoidal model.

The same annuli as those used for the radial profile were adopted for the fitting. For each annulus, the data were averaged over one beam FWHM to minimize the effects of spatial correlations. The package emcee (Foreman-Mackey et al. 2019) was used to sample the posterior distributions of $\{\delta I, \phi_{s,gas}\}$ for each annulus, using 1024 walkers taking 20,000 steps, with the first 10,000 being discarded for burn-in. Uniform priors were adopted for both parameters, with $\delta I$ assumed to be positive and less than $1 \text{ mJy beam}^{-1} \text{ km s}^{-1}$, with $0^\circ \le \phi_{s,gas} < 360^\circ$, and no annulus-to-annulus correlation. The results and their associated uncertainties were taken to be the median and 16th–84th percentile range of the marginalized posterior distributions.

The results are shown in Figure 2. Due to the substantial substructure detected in the inner disk (e.g., those discussed in

<figure><img src="image_7.png" /><figcaption>Figure 2. Results of the sub-millimeter shadow extraction. Panel (a) shows a polar deprojection of the zeroth moment map residuals shown in Figure 1(c). The error bars show the location of $\phi_{s,gas}$. The outer spiral arm reported by Teague et al. (2019b) is plotted with a dashed line with a $180^\circ$ shift showing that it aligns with the structure within $3''$. Panel (b) shows the amplitude of the variation, while panel (c) shows this as a relative fraction of the background average.</figcaption></figure>
Teague et al. (2019b), only regions beyond 1"/2 (90 au) are displayed. Inside this radius Equation (1) is a poor description of the azimuthal profile of the integrated intensity. Azimuthal variations of up to ≈0.4 mJy beam⁻¹ km s⁻¹ are found across the outer disk, relating to fractional variations of ~1% inside of 3", and increasing to up to ~30% outside 3", as the intensity rapidly drops (see Figure 1(a)). The position angle of the shadow appears to have a predominantly easterly direction (φs,gas ~ 90°), aside from a distinct deviations between 2" and 3", where φs,gas is likely dominated by the spiral features described in Teague et al. (2019b), as can be shown by the black dashed line that shows the spiral peak location (shifted by 180° to account for the fact φs,gas measures the azimuthal minimum). The reason for the spiral structure only dominating the integrated intensity out to 3" is discussed in Section 6.

## 4. Dust Shadows

HST observations have shown that a shadow is present on the surface of the disk at optical and NIR wavelengths (Debes et al. 2017). Using multi-epoch observations, Debes et al. (2017) proposed that the shadow moves in an anticlockwise direction (increasing PA), at a rate of 22°/7 yr⁻¹, equivalent to the orbital frequency of a body at 5.9 au. The authors proposed that a warped inner disk could be capable of casting such a shadow. We developed a model of the temporal evolution of the dust shadow in order to provide a comparison to the gas shadow described in Section 3.

We used the data presented in Debes et al. (2017), supplemented with an additional epoch of data from 2021 June (MJD = 59372) from the on-going HST program #16228 (J. Debes, private communication). These unpublished data were reduced following the same procedure as described in Debes et al. (2017), and all uncertainties were assumed to be 10° (J. Debes et al., 2022, in preparation). Figure 3 presents a summary of the data, with the top row showing the radial dependence of the dust shadow for the six epochs where it was spatially resolved, while the bottom panel shows the temporal evolution.

### 4.1. Curvature

It is clear from the top row of panels in Figure 3 that there is a large variation in radial morphology of the dust shadow between epochs: sometimes the shadow appears to have no or minimal radial dependence (1999, 2015, and 2016), while at other times there is a distinct radial structure (2001, 2004, and 2021). As the disk is believed to rotate in a clockwise direction (i.e., decreasing PA; Teague et al. 2019b), two of these epochs show an apparent “trailing” spiral morphology (2004 and 2021), while the 2001 data show a distinctly kinked morphology, with the outer regions hinting at a “leading” spiral morphology.

Kama et al. (2016) have proposed that light travel time effects can produce a curved (in a trailing direction) shadow if the object (called henceforth a “rim”, although the exact morphology of the obscuring feature is unknown) that casts the shadow moves an appreciable distance during the time it takes for light to reach the outer disk where the shadow is observed. For a geometrically thin, face-on disk, appropriate for TW Hya, Kama et al. (2016) showed that the radial form of the shadow is

given by

$$ \phi_{s,dust}(r, t) = \phi_{\text{rim}}(t) - \Omega_s \cdot \left( \frac{r - r_{\text{rim}}}{c} \right), \quad (2) $$

where Ωs is the orbital frequency of the shadow, c the speed of light, and (rrim, φrim(t)) describes the cylindrical coordinates of the shadowing rim. Adopting Ωs = 22°/7 yr⁻¹, the apparent orbital frequency of the shadow, we find dφs,dust/dr ≈ 4 × 10⁻⁴ deg au⁻¹, yielding variations of ≲0°/1 across the radial extent of the shadow, far too small to account for the curvature observed in the 2004 and 2021 epochs.

Alternatively, ignoring the previously constrained Ωs from the overall variation of the shadow with time, we instead assumed that Ωs is determined by the spatial morphology of the shadow on an epoch-by-epoch basis. Here, we fixed the Ωs to be the Keplerian frequency at the location of a shadowing rim, Ωs = √GM*/rrim³, and solved for rrim. Determining rrim for each of the six spatially resolved epochs yielded values of rrim ranging between rrim = 0.04 ± 0.01 au for 2004 to rrim = 0.8 ± 0.6 au for 2015. Such a vast variation in rrim is hard to reconcile with a model of a precessing inner disk (e.g., Facchini et al. 2018; Nealon et al. 2019, 2020), and suggests that the inner disk may be highly variable.

It was also argued in Kama et al. (2016) that the flared scattering surface could make a shadow appear curved, both because of the additional light travel time (the distance the light travels is the spherical polar radius of the disk, rather than the cylindrical radius), and the projection of the shadow on the sky plane. The additional distance due to a flared disk will only increase the light travel time by ≈5% for a scattering surface described by z/r = 0.3 and is thus unable to account for the observed curvature. To explore the possibility of projection effects creating curvature, we used the Python package GoFish (Teague 2019a) to model a range of flared emission surfaces and found that, given the low inclination of TW Hya, a flared scattering surface could only account for azimuthal variations of ≲0°/1 across the radial extent of the shadow.

The clear, temporally dependent, structure observed in the scattered light shadow is therefore likely to be reflective of a complex and highly variable morphology for the structure casting the shadow. As a comprehensive model of the inner disk of TW Hya is not the focus of this paper, a radially constant shadow is assumed for the rest of this analysis. Continued monitoring of the shadow at optical and NIR wavelengths will be essential to characterize the temporal variations and aid in unraveling the inner disk structure.

### 4.2. Temporal Evolution

In addition to the change in radial morphology of the shadow, the shadow as a whole is observed to move across the surface of the disk. In order to compare the gas and dust shadows, a model must be used to infer where the dust shadow lies at the time of the ALMA observations.<fnref n="10" />

First, we verify the angular frequency of the shadow under the assumption of orbital motion by fitting the temporal

<footnote n="10">¹⁰ As the ALMA observations consist of two sets of observations corresponding to the two different antenna configurations used, the date halfway between both runs is adopted as the time of the ALMA observations: 2021 July 1.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Summary of the temporal evolution of the shadow. The top row shows the six epochs where the shadow has been radially resolved in HST data (Debes et al. 2017), all sharing the same x-axis scaling. The bottom panel shows the azimuthal location of the shadow as a function of time. Two different models of the time evolution of the shadow are shown: a regular orbit, Equation (3), shown by the dashed line, or oscillatory motion, Equation (4), shown by the dotted line. The shaded regions about both lines show ±3σ about the median value of the posterior. Square points show epoch-median shadows observed in the scattered light, with each symbol representing a separate radial bin. The triangle symbol shows the median shadow location beyond 3" observed in the gas with ALMA.</figcaption></figure>

evolution of the location of the shadow with the form,

$$ \tilde{\phi}_{s,dust}(t) = (\bar{\Omega}_s \cdot (t - \delta t)) \mod 360^{\circ} \quad (3) $$

where $\bar{\Omega}_s$ is the angular frequency of the shadow and $\delta t$ is a phase offset. As in Section 3, *emcee* was used to sample the posterior distributions of {$\bar{\Omega}_s, \delta t$} using 1024 walkers taking 20,000 steps, with the first 10,000 being discarded for burn-in. The 50th percentile and 16th–84th percentile range of the posterior distributions were adopted as the best-fit value and uncertainty, respectively. We stress that these uncertainties only represent the statistical uncertainty on these parameters, and do not reflect the systematic uncertainties associated with model selection. The angular frequency of the shadow was found to be $\bar{\Omega}_s = 22^{\circ}9 \pm 0.1 \text{ yr}^{-1}$, consistent with the angular frequency reported in Debes et al. (2017), and shown in the bottom panel of Figure 3 as a dashed line.

Both the 2004 data and the new epoch (2021) of HST data appear to show a shadow that is at a smaller PA than would be predicted by this orbital motion, as marked by the dashed lines in Figure 3. This motivates exploration of an alternative model of the temporal evolution of the shadow. One possibility is oscillatory motion, due perhaps to a differentially precessing inner disk (e.g., Facchini et al. 2018; Nealon et al. 2019, 2020). To check the suitability of this model, we fit the data with the

form

$$ \tilde{\phi}_{s,dust}(t) = \delta\tilde{\phi}_s \sin\left(2\pi \frac{t - \delta t}{\tau}\right) + \langle\tilde{\phi}_{s,dust}\rangle, \quad (4) $$

where $\langle\tilde{\phi}_{s,dust}\rangle$ is the average position angle of the shadow, $\delta\tilde{\phi}_{s,dust}$ is the amplitude of the oscillation, $\tau$ is the period of the oscillation, and $\delta t$ is a phase offset. Note that a bar is used to denote the assumption of orbital motion (constant angular velocity of the shadow), as in Equation (3), while a tilde represents the oscillatory motion of the shadow, as described by Equation (4).

Using the same fitting procedure as before, we found $\delta\tilde{\phi}_{s,dust} = 55^{\circ} \pm 2^{\circ}$, $\langle\tilde{\phi}_{s,dust}\rangle = 239^{\circ} \pm 2^{\circ}$, and $\tau = 15.8 \pm 0.1 \text{ yr}$. This fit is shown by the dotted line in the bottom panel of Figure 3. To verify that this fit is not biased by the additional epoch of observations, we fit only those epochs reported in Debes et al. (2017) and find values of $\delta\tilde{\phi}_{s,dust} = 57^{\circ} \pm 2^{\circ}$, $\langle\tilde{\phi}_{s,dust}\rangle = 232^{\circ} \pm 3^{\circ}$, and $\tau = 15.8 \pm 0.1 \text{ yr}$. The similarity in inferred parameters suggests that the original data sets already support the oscillatory model.

Unsurprisingly, both models give a similar period of ~15.8 yr due to the timing of the various observations. However, while the oscillatory motion is unable to describe the 2005 data, potentially owing to the lower signal-to-noise of the data compared to the other epochs (J. Debes, private
communication), it appears to give a better fit to both the 2004 and 2021 epochs. More specifically, for 2004 the shadow was observed at $\phi_{s,dust} = 288^{\circ} \pm 37^{\circ}$ (where the quoted uncertainty describes the uncertainty-weighted standard deviation of the sample), while the predicted locations were $\tilde{\phi}_{s,dust} = 319^{\circ} \pm 1^{\circ}$ and $\tilde{\phi}_{s,dust} = 288^{\circ} \pm 2^{\circ}$, for the orbital and oscillatory motion, respectively, with the oscillatory motion finding a much better agreement. In a similar manner the 2021 epoch data showed the shadow located at $\phi_{s,dust} = 298^{\circ} \pm 20^{\circ}$, for which the orbital model predicts a location of $\tilde{\phi}_{s,dust} = 350^{\circ} \pm 2^{\circ}$ while the oscillatory model predicts $\tilde{\phi}_{s,dust} = 281^{\circ} \pm 3^{\circ}$, again resulting in a better fit for the oscillatory motion. While these data are suggestive that oscillatory motion would better describe the temporal evolution of the shadow, additional epochs are required to unambiguously distinguish between these scenarios and, as such, both models are used in the remainder of this work.

## 5. Offsets between Gas and Dust Shadows

Although the mass of a protoplanetary disk is dominated by the gas, it is the dust that determines the temperature. Dust grains absorb the stellar radiation which then transfers the heat to the surrounding gas through collisions between gas molecules and dust grains. When a disk region enters a shadow, there is a deficit of photons for the dust grains to absorb and to scatter, resulting in a drop in the dust temperature. The gas may not immediately experience a drop in temperature: gas molecules must first collide with the cooler dust grains to transfer some of their excess energy and cool themselves, described by the thermal accommodation timescale, $t_{th}$. This timescale combines the collisional rate of gas and dust molecules, $t_{coll}$, and the speed at which the collisional energy can be converted to internal thermal energy, and is given by

$$t_{th} = \frac{c_V}{2k_B\tilde{\alpha}_T} t_{coll}, \qquad (5)$$

where $c_V$ is the specific heat capacity of the gas, $k_B$ is the Boltzmann constant, $\tilde{\alpha}_T \approx 0.5$ and is the thermal accommodation coefficient that characterizes the efficiency of the heat transfer between gas molecules and dust grains (Burke & Hollenbach 1983). The prefactor is sufficiently close to unity that $t_{th} \approx t_{coll}$, is assumed implicitly for the remainder of this paper.

Indeed, it has been shown that the gas-dust collisional timescale in the surface layers of protoplanetary disks may not be infinitely short due to low dust abundance (Facchini et al. 2017; Bae et al. 2021). If $t_{coll}$ is a reasonable fraction of the local orbital period, $t_{orb}$, then an offset is expected between the shadowed location traced by scattered light and the region where the gas is coldest as the gas takes time to react to the change in dust temperature.

Using the models of the temporal evolution of the dust shadows described in the previous section, we calculated the azimuthal location of the shadow at the time of the ALMA observations. As Debes et al. (2017) does not discuss the presence or absence of shadows beyond $2''/4$, we have assumed that the shadow persists to the outer disk where the gas shadow is detected, and modeled the azimuthal offset between dust and gas at different radii assuming both the orbital and oscillatory models. We found an azimuthal offset between the gas and dust

<figure><img src="image_9.png" /><figcaption>Figure 4. Comparing the gas shadow, black points, with the predicted dust shadow location for both oscillatory and orbital motion, blue and red points, respectively. For each dust shadow model both a straight shadow and curved shadow model are shown with circle and square markers, respectively. The background image is a de-saturated version of Figure 1(c). Note that disk is rotating in a clockwise direction.</figcaption></figure>

shadows of $\sim 60^{\circ}-0^{\circ}$ between $3''/1$ and $3''/9$, as shown in Figure 4. To verify that the gas and dust shadows cannot be comoving, we repeated the fitting described in Section 4.2 including the ALMA data point. There was negligible change in the resulting model parameters, confirming that neither the orbital nor the oscillatory model could simultaneously explain the temporal evolution of the gas and dust shadows and that an offset is required, except for the outermost part of the disk. The radially resolved angular separation between shadows is shown in Figures 5(a) and (c) for a straight and curved shadow, respectively. In each panel both the orbital and oscillatory models of the shadows' movement are shown in red and blue, respectively. We propose that this offset could be due to the finite gas-dust collisional timescale.

The angular separation between the shadows in the dust and the gas can therefore be written as

$$\phi_{s,dust}(r) - \phi_{s,gas}(r) = (\Omega_s - \Omega_{\text{kep}}(r)) \cdot t_{\text{coll}}, \qquad (6)$$

where $\Omega_s - \Omega_{\text{kep}}(r)$ describes the relative angular velocity between the shadow and a parcel of gas in the disk rotating at the Keplerian frequency. Note that while $\Omega_{\text{kep}}$ varies radially, $\Omega_s$ is constant as a function of radius as it only depends on the angular velocity of the shadowing rim. At the time of the ALMA observations, both models of the shadows' movement predict that the shadow is moving counterclockwise on the sky (increasing $\phi$), in the opposite direction to the rotation of the disk such that $\partial\phi_{s,dust}/\partial t$ is positive and $\Omega_{\text{kep}}$ is negative. As $\Omega_s$ is considerably larger than $\Omega_{\text{kep}}$, any changes in the gas velocity structure, such as substantial sub-Keplerian rotation arising due to a large, negative pressure gradient at the disk edge (e.g., Dullemond et al. 2020), would have negligible impact on this result.
<figure><img src="image_10.png" /><figcaption>Figure 5. Angular separation between dust and gas shadows and the resulting inferred $t_{coll}$ assuming a straight and curved dust shadow, left and right columns, respectively. In each column, the top panels, (a) and (c), show the angular separation between the predicted location of the shadow traced by scattered light and traced by CO emission in the outer disk. Red points assume an orbital motion, while blue points assume oscillatory motion. In each column, the bottom panels, (b) and (d), show the calculated $t_{coll}$ for both temporal models, with dashed lines in the background showing contours of constant $t_{coll}/t_{orb}$. In both panels the error bars represent $1\sigma$ uncertainties, while arrows represent upper limits. A small radial offset has been applied to both points to avoid overlapping error bars.</figcaption></figure>

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

Adopting the same stellar mass as used for the Keplerian CLEAN mask in Section 2 and the instantaneous angular velocity of the shadow at the time of observations, $d\tilde{\phi}_{s,dust}/dt = 22.6 \text{ yr}^{-1} \pm 0.1 \text{ yr}^{-1}$ and $d\tilde{\phi}_{s,dust}/dt = 5.1 \text{ yr}^{-1} \pm 0.8 \text{ yr}^{-1}$, the angular separation of shadows can be related to $t_{coll}$, shown in Figure 5(c). Dashed lines in the background of this panel show lines of constant $t_{coll}/t_{orb}$ ($t_{coll}$ is given by Equation (7)), where values range between roughly $10^{-4}$ and $10^{-2}$. Beyond $3''5$, the posteriors of $t_{coll}$ are consistent with instantaneous gas-dust collisions for the oscillatory model. We note that it is the much slower angular velocity of the shadow if oscillatory motion is assumed than for the orbital motion that gives rise to a larger $t_{coll}$, despite a smaller angular separation between the gas and dust shadows.

There appears to be a subtle radial decrease in $t_{coll}$ for both models, suggesting potentially a more efficient coupling between gas and grains in the outer disk for the vertical layer probed by CO emission, an idea discussed further in Section 6.1. An alternative scenario is that some of the radial dependency arises from curvature in the dust shadow (see Section 4.1) that was not accounted for in the angular offset between gas and dust shadows. To understand the impact of a curved dust shadow, we repeated the calculation of $t_{coll}$ under the assumption of a curved dust shadow. We adopted a radial morphology given by Equation (2) (i.e., an arithmetic spiral), and linearly interpolated between the best-fit models of the 2016 and 2021 epochs for $\phi_{s,dust}$ to model the curvature of the shadow during the ALMA observations. This results in a morphology with a radial dependence of $d\tilde{\phi}_{s,dust}/dr \approx 0.29 \text{ au}^{-1}$, such that average offset between shadows are increased by approximately $70^\circ$, as shown in Figure 5(c). A large angular offset between shadows results in an increase of $t_{coll}$ to ~5 yr for the orbital motion and ~15 yr for the oscillatory motion (see Figure 5(d)). A slight radial dependence

of $t_{coll}$ potentially persists; however, these results suggest that a curved dust shadow could equally well explain any radial trends.

## 6. Discussion

We have shown that CO $J=2-1$ emission from the disk around TW Hya shows azimuthal variations on the order of 20% in the outer disk ($r \gtrsim 3''$), likely related to a shadow observed in scattered light at smaller radii ($r \lesssim 2''5$; Debes et al. 2017). Modeling the temporal evolution of the dust shadow suggests an offset between the gas and dust shadow which was interpreted as due to a non-negligible collisional timescale, $t_{coll}$. In this section we discuss the implications of this finding.

### 6.1. Collisional Timescale

Following Bae et al. (2021), the collisional timescale for gas molecules is given by

$$ 
\begin{aligned}
t_{\text{coll}} &= \left(\frac{4}{3}\rho_s\right) \left(\frac{\langle a^3 \rangle}{\langle a^2 \rangle}\right) \left(\frac{m_g}{\rho_g v_{\text{th}}}\right) \left(\frac{\rho_d}{\rho_g}\right)^{-1}, \\
&\approx 30.4 \text{ yr} \left(\frac{\rho_s}{3 \text{ g cm}^{-3}}\right) \left(\frac{\langle a^3 \rangle / \langle a^2 \rangle}{0.1 \text{ µm}}\right) \left(\frac{T_g}{20 \text{ K}}\right)^{-1/2} \\
&\quad \times \left(\frac{\rho_g}{10^{-15} \text{ g cm}^{-3}}\right)^{-1} \left(\frac{\rho_d / \rho_g}{0.01}\right)^{-1},
\end{aligned} 
\quad (7) $$

where $\rho_s$ is the bulk density of the dust grains, $\langle a^3 \rangle / \langle a^2 \rangle$ is the mean dust grain size, $m_g$ is the mean mass of gas molecules, $\rho_g$ is the local gas volume density, $v_{th}$ is the thermal velocity of the gas, and $\rho_d / \rho_g$ is the local dust-to-gas ratio. We note that $t_{coll}$ therefore scales with the local mean grain size, the local gas-to-
dust ratio, and inversely with the square root of the gas temperature (through the thermal velocity of the gas; see also Young et al. 2004; Facchini et al. 2017).

With some representative numbers assuming that the line is marginally optically thin, such that it arises from colder regions closer to the midplane, i.e., $\rho_s = 3 \text{ g cm}^{-3}$, $\langle a^3 \rangle / \langle a^2 \rangle = 0.1 \mu\text{m}$, $\rho_g = 10^{-15} \text{ g cm}^{-3}$ (an estimate for one scale height above the midplane assuming a surface density $\Sigma_g = 1 \text{ g cm}^{-2}$ and a gas scale height of 20 au, broadly consistent with the model of Calahan et al. 2021), $T_g = 20 \text{ K}$, and $\rho_d/\rho_g = 0.01$, the collisional timescale is $t_{\text{coll}} \simeq 30.4 \text{ yr}$. This representative value is several times larger than the $t_{\text{coll}}$ inferred from the observed shadow offset in TW Hya (bottom panels of Figure 5). Such a discrepancy could easily be accounted for with variations of a factor of a few for all the properties described in Equation (7). In particular, $\rho_d/\rho_g$ can be one to a few orders of magnitude smaller than 0.01 in surface layers when the settling of dust grains is significant (see, e.g., Figure 2 of Bae et al. 2021), as would be the case for the low levels of turbulence inferred for this disk (Hughes et al. 2011; Teague et al. 2016; Flaherty et al. 2018).

## 6.2. Excitation Conditions

The shadows described in this work were not observed in previously published high-resolution ALMA data of CO emission (Huang et al. 2018; Teague et al. 2019b). Instead, the previous observations displayed spiral arms that were observed to extend to larger radii (just outside 200 au or $3''$4), before the CO emission was no longer detectable. The difference between these archival observations and those presented here are that the former were at slightly higher frequencies, observing the $J=3-2$ transition with $E_u = 33.2 \text{ K}$, as opposed to the lower energy $E_u = 16.6 \text{ K}$ of the $J=2-1$ emission presented here. It is therefore probable that the excitation conditions of the emission influence where and when a shadow can be observed.

In general lower-energy transitions will probe deeper layers in the disk, tracing regions closer to the midplane, than higher-energy transitions (e.g., Dartois et al. 2003). As the vertical settling of dust grains increases with height, such that the local dust-to-gas ratio drops, higher-frequency transitions will probe a region of the disk with longer $t_{\text{coll}}$. When $t_{\text{coll}}$ increases to become a substantial fraction of $t_{\text{orb}}$ then a gas parcel will move out of the shadowed region before it has time to collide with the dust to cool down and create a shadow. Therefore, if the $J=3-2$ emission were sampling a more elevated region in the disk having a longer $t_{\text{coll}}$ than the shadow-crossing time this could explain why shadows were not seen, but spirals were (which Bae et al. 2021 and Muley et al. 2021 argue should get stronger with elevation).

To explore this scenario, we need to consider when a shadow is observable in molecular line emission. The condition for a shadow to form is given by

$$t_{\text{coll}}(r) < \frac{\Delta\phi_{s,\text{dust}}(r)}{|\Omega_s - \Omega_{\text{kep}}(r)|}, \quad (8)$$

where $\Delta\phi_{s,\text{dust}}(r)$ is the radially dependent azimuthal extent of the dust shadow, such that the right-hand side is the shadow crossing time. We adopt a radially constant $\Delta\phi_s \simeq 180^\circ$ for the shadow and the maximum angular velocity of the shadow, $\Omega_s = 22.9 \text{ yr}^{-1}$ (as any smaller $\Omega_s$ will just increase the limits

of $t_{\text{coll}}$). As $|\Omega_s| \gg |\Omega_{\text{kep}}|$ at these radii, it is irrelevant if the disk is rotating clockwise or counterclockwise, as the relative motion is dominated by the movement of the shadow. We find $t_{\text{coll}} \lesssim 8 \text{ yr}$ is required for a gas to feel the effects of the shadow. A shadow crossing time of $\sim 8 \text{ yr}$ sits squarely in the range of $t_{\text{coll}}$ inferred for CO $J=2-1$ considering the different radial and temporal dependencies of the dust shadow. Thus, it is likely the $J=2-1$ emission probes a vertical layer where $t_{\text{coll}}$ is close to this limit, and $J=3-2$ probes an elevated layer where $t_{\text{coll}} > 8 \text{ yr}$ due to the decrease in gas and dust density.

This scenario also provides an explanation for why spirals are clearly seen with the $J=2-1$ emission but are quickly lost beyond $r \sim 3''$. As the surface density drops toward the outer disk, the optical depth of the line decreases such that the emission probes a progressively less elevated region (which has been directly observed in disks before, e.g., Teague et al. 2019a; Law et al. 2021). This would result in the emission tracing progressively smaller $t_{\text{coll}}$ values until a threshold is reached (Equation (8)) and shadows can be observed in the gas. A direct test of this scenario would be to search for similar shadows or azimuthal variations in emission which arise from deeper regions of the disk, either the $J=1-0$ transition of CO, or emission from less abundant isotopologues, and to identify where the shadows first become observable, although currently no observations of such emission at a sufficient sensitivity or angular resolution exist.

## 7. Summary

In this paper we present new observations of CO $J=2-1$ emission from the disk around TW Hya. The integrated intensity exhibits azimuthal variations in the outer disk reaching fractional variations of up to $\sim 20\%$ in an east/west direction. These variations are consistent in structure with the azimuthal asymmetries and shadows detected previously in scattered light observations (Roberge et al. 2005; Debes et al. 2017).

To ascertain any offset in shadows in the dust and gas, observed with HST and ALMA, respectively, the temporal evolution of the shadow was re-analyzed. It was found that an oscillatory motion of the dust shadow, rather than orbital motion, would better explain the temporal evolution and readily account for the apparent counter-rotation movement of the shadow over the face of the disk. Both models of the movement of the shadow resulted in a $10^\circ-60^\circ$ offset between the dust and gas shadows which was interpreted as due to a non-negligible collisional timescale between the gas molecules and dust grains. Given the uncertainty in the angular velocity of the shadow, this corresponds to a $t_{\text{coll}}$ of between 1 and 10 yr.

Although the HST data are suggestive of a curved shadow, it is unclear with the current data what could cause this curvature —both light travel time and projection effects proposed by Kama et al. (2016) can be ruled out. Adopting a simple linear model of the shadow curvature in the analysis demonstrated that a curved dust shadow would increase the angular separation between dust and gas shadows, resulting in a lengthened $t_{\text{coll}}$, ranging between 5 and 20 yr.

The inferred $t_{\text{coll}}$ timescales are close to the shadow crossing time, suggesting that a small change in local physical conditions would be sufficient for the gas to no longer exhibit a shadow. This is likely the reason the gas shadow is not seen in previous observations of $J=3-2$ emission as the emission
originates from a more elevated region where $t_{\text{coll}}$ is larger due to the lower dust densities.

Software: emcee (Foreman-Mackey et al. 2019), GoFish (Teague 2019a), bettermoments (Teague & Foreman-Mackey 2018), CASA (v5.8.0; McMullin et al. 2007), keplerian_mask.py (Teague 2020).

The authors would like to thank the referee, Takayuki Muto, for their careful and helpful review of the manuscript, John Debes for graciously sharing unpublished HST data from their program #16228, and enlightening discussion with Kees Dullemond. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.A.00021.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with program #16228. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 101002188) Support for J.H. was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51460.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

## ORCID iDs

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

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Stefano Facchini https://orcid.org/0000-0003-4689-2684

Jane Huang https://orcid.org/0000-0001-6947-6072

David Wilner https://orcid.org/0000-0003-1526-7587

## References

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

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

* Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

* Apai, D., Pascucci, I., Brandner, W., et al. 2004, A&A, 415, 671

* Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56

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

* Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223

* Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8

* Canta, A., Teague, R., Le Gal, R., & Öberg, K. I. 2021, ApJ, 922, 62

* Casassus, S., Pérez, S., Osses, A., & Marino, S. 2019, MNRAS, 486, L58

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

* Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2

* Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773

* Debes, J. H., Jang-Condell, H., & Schneider, G. 2016, ApJL, 819, L1

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

* Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

* Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137

* Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16

* Facchini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459

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

* Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, JOSS, 4, 1864

* Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

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

* 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

* Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

* Kama, M., Pinilla, P., & Heays, A. N. 2016, A&A, 593, L20

* Krist, J. E., Stapelfeldt, K. R., Ménard, F., Padgett, D. L., & Burrows, C. J. 2000, ApJ, 538, 793

* Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4

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

* Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33

* McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

* Muley, D., Dong, R., & Fung, J. 2021, AJ, 162, 129

* Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951

* Nealon, R., Price, D. J., & Pinte, C. 2020, MNRAS, 493, L143

* 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

* Roberge, A., Weinberger, A. J., & Malumuth, E. M. 2005, ApJ, 622, 1171

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

* Teague, R. 2019a, JOSS, 4, 1632

* Teague, R. 2019b, JOSS, 4, 1220

* Teague, R. 2020, richteague/keplerian_mask: Initial Release, v1.0, Zenodo doi: 10.5281/zenodo.4321137

* Teague, R., Bae, J., & Bergin, E. A. 2019a, Natur, 574, 378

* Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019b, ApJL, 884, L56

* Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

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

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

* Walsh, C., Loomis, R. A., Öberg, K. I., et al. 2016, ApJL, 823, L10

* Weinberger, A. J., Becklin, E. E., Schneider, G., et al. 2002, ApJ, 566, 409

* Young, K. E., Lee, J.-E., Evans, N. J. I., Goldsmith, P. F., & Doty, S. D. 2004, ApJ, 614, 252

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