# Gas and Dust Shadows in the TW Hydrae Disk

Richard Teague Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Jaehan Bae Department of Astronomy, University of Florida, Gainesville, FL 32611, USA  Myriam Benisty Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France  Sean M. Andrews Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Stefano Facchini Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy  Jane Huang NASA Hubble Fellowship Program Sagan Fellow  Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, United States of America  David Wilner Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

###### Abstract

We present new observations of CO $J=2-1$ emission from the protoplanetary
disk around TW Hya. Emission is detected out to $240~{}{\rm au}$ ($4\arcsec$)
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
reanalysis of the multi-epoch observations of the dust shadow in scattered
light from Debes et al. (2017) 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 ${\sim}~{}100\arcdeg$. 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_{\rm coll}$. The angular offsets derived are equivalent to collisional
timescales that range from the near instantaneous up to $t_{\rm coll}\sim
10~{}{\rm years}$, depending on whether a straight or curved dust shadow, as
suggested by HST observations reported by Debes et al. (2017), is adopted. The
inferred range of $t_{\rm 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_{\rm coll}$.

††journal: ApJ††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)

##  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\approx 5\fdg 8$;
Teague et al., 2019b). As such, it is often used as a young analogue 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-) mm wavelengths, tracing the thermal emission of the mm-
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-µm 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 au 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-mm 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 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\fdg 7~{}{\rm yr^{-1}}$ in an anti-clockwise 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\farcs 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 reanalysis 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 12CO $J=2-1$,
CN $N=2-1$ and CS $J=5-4$. The correlator was set up to cover the 12CO and CS
line with a channel spacing of 15.3 kHz (resulting in a 30.5 kHz resolution,
equivalent to $40~{}{\rm m\,s^{-1}}$ and 230.538 GHz), while the CN was
observed at a coarser 61 kHz channel spacing (yielding a 122 kHz resolution,
equivalent to $160~{}{\rm m\,s^{-1}}$) 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 CH2CN in the continuum window was presented in Canta et al.
(2021). In this paper we focus on the CO emission, leaving analysis of the CS
and CN for future work.

###  2.1 Calibration

Figure 1: Summary of the CO integrated intensity. (a) The 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) The
$J=2-1$ moment map with the beam size shown by the ellipse in the lower left
corner. The solid contours show SNR of $5\times 2^{n}$ with
$n=\\{0,\,1,\,\dots,\,7\\}$. (c) The 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 and $5~{}\sigma$ contours. The dashed
line is a guide to emphasize the global east/west asymmetry.

Short spacing data (baselines of 15 m to 500 m) were taken 2019 April 4th with
two executions including 47.1 minutes of on-source time each. During these
executions the quasars J1037-2934 and J1147-3812 were observed for calibration
purposes, the former acted 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 Sep 29th. Only 4 out of the requested 12 executions
were observed due to scheduling constraints. Each execution included 50.1
minutes 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 hours.

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~{}{\rm m\,s^{-1}}$ (the spectral
resolution of the data) was 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\farcs 19\times 0\farcs 17$ ($89\arcdeg$). A Keplerian
mask was generated,111github.com/richteague/keplerian_mask making sure that
all disk emission was contained within the mask. Following Teague et al.
(2019b), a stellar mass of $M_{\star}=0.81~{}M_{\odot}$ and a viewing geometry
described by $i=5\fdg 8$ and ${\rm PA}=151\arcdeg$ 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~{}{\rm mJy~{}beam^{-1}}$. Integrating over the Keplerian
mask we recover an integrated intensities of $18.31\pm 0.03~{}{\rm
Jy~{}km\,s^{-1}}$. This uncertainty does not include a systematic uncertainty
of ${\sim}~{}10\%$ associated with the flux calibration of the data.

The package bettermoments222github.com/richteague/bettermoments was used to
calculate the zeroth moment map of the data. This was calculated using the
Keplerian CLEAN mask without applying any $\sigma$-clipping. The resulting
moment map is shown in Fig. 1b. The analysis of additional moment maps is the
focus of a companion paper (Teague et al., in prep.).

##  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: $\Delta{\rm RA}=20.0\pm 0.1~{}{\rm
mas}$ and $\Delta{\rm Dec}=2.0\pm 0.1~{}{\rm 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 was binned into annuli with a width of 1/4
of the beam FWHM (${\approx}~{}50~{}{\rm 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 Fig.
1a, 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 Fig 1b.
The residuals are shown in Fig. 1c. 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\arcsec$, 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 misspecified center, the same
residuals were calculated for a range of offsets 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 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\big{(}\phi-\phi_{s,\,{\rm
gas}}(r)-\pi\big{)}+\langle I(r)\rangle,$ (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,\,{\rm 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,\,{\rm 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,\,{\rm gas}}$ does not measure the azimuthal
minimum of each annulus, but rather than 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 was 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,\,{\rm 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~{}{\rm mJy\,beam^{-1}~{}km\,s^{-1}}$, while
$0\arcdeg\leq\phi_{s,\,{\rm gas}}<360\arcdeg$, and no annulus-to-annulus
correlation. The results and their associated uncertainties were taken to be
the median and 16th to 84th percentile range of the marginalized posterior
distributions.

The results are shown in Fig. 2. Due to the substantial substructure detected
in the inner disk (e.g., those discussed in Teague et al., 2019b), only
regions beyond $1\farcs 2$ (90 au) are displayed. Inside this radius Eqn. 1 is
a poor description of the azimuthal profile of the integrated intensity.
Azimuthal variations of up to ${\approx}~{}0.4~{}{\rm
mJy\,beam^{-1}~{}km\,s^{-1}}$ are found across the outer disk, relating to
fractional variations of ${\sim}~{}1\%$ inside of 3″, and increasing to up to
${\sim}~{}30\%$ outside 3″, as the intensity rapidly drops (see Fig. 1a). The
position angle of the shadow appears to have a predominantly easterly
direction ($\phi_{s,\,{\rm gas}}\sim 90\arcdeg$), aside from a distinct
deviations between 2″ and 3″, where $\phi_{s,\,{\rm 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 $\phi_{s,\,{\rm gas}}$ measures the azimuthal minimum).
The reason for the spiral structure only dominating the integrated intensity
out to 3″ is discussed in Section 6.

Figure 2: Results of the sub-mm shadow extraction. Panel (a) shows a polar
deprojection of the zeroth moment map residuals shown in Fig. 1c. The error
bars show the location of $\phi_{s,\,{\rm gas}}$. The outer spiral arm
reported by Teague et al. (2019b) is plotted with a dashed line with a 180°
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.

##  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 anti-
clockwise direction (increasing position angle), at a rate of $22\fdg 7~{}{\rm
yr^{-1}}$, 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 June 2021 (MJD = 59372) from the on-going HST
program #16228 (J. Debes, private communication). This unpublished data was
reduced following the same procedure as described in Debes et al. (2017), and
all uncertainties were assumed to be $10\arcdeg$ (Debes et al., in prep.).
Figure 3 presents a summary of the data, with the top row showing the radial
dependence of the dust shadow for the 6 epochs where it was spatially
resolved, while the bottom panel shows the temporal evolution.

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, Eqn. 3, shown
by the dashed line, or oscillatory motion, Eqn. 4, shown by the dotted line.
The shaded regions about both lines show $\pm 3\sigma$ 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.

###  4.1 Curvature

It is clear from the top row of panels in Fig. 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 shows 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 if 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,\,{\rm dust}}(r,t)=\phi_{\rm
rim}(t)-\Omega_{s}\cdot\left(\frac{r-r_{\rm rim}}{c}\right),$ (2)

where $\Omega_{s}$ is the orbital frequency of the shadow, $c$ the speed of
light, and $(r_{\rm rim},\,\phi_{\rm rim}(t))$ describing the cylindrical
coordinates of the shadowing rim. Adopting $\Omega_{s}=22\fdg 7~{}{\rm
yr^{-1}}$, the apparent orbital frequency of the shadow, we find ${\rm
d}\phi_{s,\,{\rm dust}}/{\rm d}r\approx 4\times 10^{-4}~{}{\rm
deg~{}au^{-1}}$, yielding variations of $\lesssim 0\fdg 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 $\Omega_{s}$ from the
overall variation of the shadow with time, we instead assumed that
$\Omega_{s}$ is determined by the spatial morphology of the shadow on an
epoch-by-epoch basis. Here, we fixed the $\Omega_{s}$ to be the Keplerian
frequency at the location of a shadowing rim, $\Omega_{s}=\sqrt{GM_{*}/r_{\rm
rim}^{3}}$, and solved for $r_{\rm rim}$. Determining $r_{\rm rim}$ for each
of the six spatially resolved epochs yielded values of $r_{\rm rim}$ ranging
between $r_{\rm rim}=0.04\pm 0.01$ au for 2004 to $r_{\rm rim}=0.8\pm 0.6$ au
for 2015. Such a vast variation in $r_{\rm rim}$ 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 $\approx 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 $\lesssim 0\fdg 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.333As 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: July 1st, 2021.

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

$\bar{\phi}_{s,\,{\rm dust}}(t)=\big{(}\bar{\Omega}_{s}\cdot(t-\delta
t)\big{)}\,\,{\rm mod}\,\,360\arcdeg$ (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 to 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\fdg 9\pm 0\fdg 1~{}{\rm yr^{-1}}$, consistent with the
angular frequency reported in Debes et al. (2017), and shown in the bottom
panel of Fig. 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 Fig. 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,\,{\rm
dust}}(t)=\delta\tilde{\phi}_{s}\sin\left(2\pi\,\frac{t-\delta
t}{\tau}\right)+\langle\tilde{\phi}_{s,\,{\rm dust}}\rangle,$ (4)

where $\langle\tilde{\phi}_{s,\,{\rm dust}}\rangle$ is the average position
angle of the shadow, $\delta\tilde{\phi}_{s,\,{\rm 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 symbol is used to denote the assumption of
orbital motion (constant angular velocity of the shadow), as in Eqn. 3, while
a tilde represents the oscillatory motion of the shadow, as described by Eqn.
4.

Using the same fitting procedure as before, we found
$\delta\tilde{\phi}_{s,\,{\rm dust}}=55\arcdeg\pm 2\arcdeg$,
$\langle\tilde{\phi}_{s,\,{\rm dust}}\rangle=239\arcdeg\pm 2\arcdeg$ and
$\tau=15.8\pm 0.1~{}{\rm yr}$. This fit is shown by the dotted line in the
bottom panel of Fig. 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,\,{\rm
dust}}=57\arcdeg\pm 2\arcdeg$, $\langle\tilde{\phi}_{s,\,{\rm
dust}}\rangle=232\arcdeg\pm 3\arcdeg$ and $\tau=15.8\pm 0.1~{}{\rm yr}$. The
similarity in inferred parameters suggests that the original datasets already
support the oscillatory model.

Unsurprisingly, both models give a similar period of ${\sim}~{}15.8~{}{\rm
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,\,{\rm dust}}=288\arcdeg\pm 37\arcdeg$ (where the quoted uncertainty
describes the uncertainty-weighted standard deviation of sample), while the
predicted locations were $\bar{\phi}_{s,\,{\rm dust}}=319\arcdeg\pm 1\arcdeg$
and $\tilde{\phi}_{s,\,{\rm dust}}=288\arcdeg\pm 2\arcdeg$, 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,\,{\rm dust}}=298\arcdeg\pm 20\arcdeg$, for which
the orbital model predicts a location of $\bar{\phi}_{s,\,{\rm
dust}}=350\arcdeg\pm 2\arcdeg$ while the oscillatory model predicts
$\tilde{\phi}_{s,\,{\rm dust}}=281\arcdeg\pm 3\arcdeg$, 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_{\rm th}$.
This timescale combines the collisional rate of gas and dust molecules,
$t_{\rm coll}$, and the speed at which the collisional energy can be converted
to internal thermal energy, and is given by

$t_{\rm th}=\frac{c_{V}}{2k_{B}\tilde{\alpha}_{T}}~{}t_{\rm coll},$ (5)

where $c_{V}$ is the specific heat capacity of the gas, $k_{\rm B}$ is the
Boltzmann constant, $\tilde{\alpha}_{\rm 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_{\rm th}\approx t_{\rm
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_{\rm coll}$ is
a reasonable fraction of the local orbital period, $t_{\rm orb}$, then an
offset is expected between the shadowed location traced by scattered light and
the region where the gas is the coldest as the gas takes time to react to the
change in dust temperature.

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 Fig. 1c. Note that disk is
rotating in a clockwise direction.

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\farcs 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 shadows of $\sim 60\arcdeg-0\arcdeg$ between $3\farcs 1$ and
$3\farcs 9$, as shown in Fig. 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 Figs. 5a and
c for a straight and curved shadow, respectively. In each panel both the
orbital and oscillatory models of the shadows movement in red and blue,
respectively. We propose that this offset could be due to finite gas-dust
collisional timescale.

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

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

where $\Omega_{s}-\Omega_{\rm 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_{\rm 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 counter-clockwise on the
sky (increasing $\phi$), in the opposite direction to the rotation of the disk
such that $\partial\phi_{s,\,{\rm dust}}/\partial t$ is positive and
$\Omega_{\rm kep}$ is negative. As $\Omega_{s}$ is considerably larger than
$\Omega_{\rm 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 5: The angular separation between dust and gas shadows and the
resulting inferred $t_{\rm coll}$ assuming a straight and curved dust shadow,
left and right columns, respectively. In each column, the top panel, (a) and
(c), shows 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 panel, (b) and (d), shows the calculated $t_{\rm
coll}$ for both temporal models, with dashed lines in the background showing
contours of constant $t_{\rm coll}\,/\,t_{\rm 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.

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, ${\rm d}\bar{\phi}_{s,\,{\rm dust}}/{\rm d}t=22\fdg 6~{}{\rm
yr^{-1}}\pm 0\fdg 1~{}{\rm yr^{-1}}$ and ${\rm d}\tilde{\phi}_{s,\,{\rm
dust}}/{\rm d}t=5\fdg 1~{}{\rm yr^{-1}}\pm 0\fdg 8~{}{\rm yr^{-1}}$, the
angular separation of shadows can be related to $t_{\rm coll}$, shown in Fig.
5c. Dashed lines in the background of this panel show lines of constant
$t_{\rm coll}\,/\,t_{\rm orb}$ ($t_{\rm coll}$ is given by Eqn. 7), where
values range between roughly $10^{-4}$ and $10^{-2}$. Beyond $3\farcs 5$, the
posteriors of $t_{\rm 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_{\rm coll}$, despite a smaller
angular separation between the gas and dust shadows.

There appears to be a subtle radial decrease in $t_{\rm 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_{\rm coll}$ under the assumption of a curved dust shadow. We adopted a
radial morphology given by Eqn. 2 (i.e., an arithmetic spiral), and linearly
interpolated between the best-fit models of the 2016 and 2021 epochs for
$\phi_{s,\,{\rm dust}}$ to model the curvature of the shadow during the ALMA
observations. This results in a morphology with a radial dependence of ${\rm
d}\phi_{s,\,{\rm dust}}/{\rm d}{r}\approx 0\fdg 29~{}{\rm au^{-1}}$, such that
average offset between shadows are increased by approximately $70\arcdeg$, as
shown in Fig. 5c. A large angular offset between shadows results in an
increase of $t_{\rm coll}$ to ${\sim}~{}5$ yr for the orbital motion and
${\sim}~{}15$ yr for the oscillatory motion (see Fig. 5d). A slight radial
dependence of $t_{\rm 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\arcsec$), likely related to a shadow observed in scattered light at smaller
radii ($r\lesssim 2\farcs 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_{\rm 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

$\displaystyle t_{\rm coll}$ $\displaystyle=$ $\displaystyle\left({4\over
3}\rho_{s}\right)\left({\langle a^{3}\rangle\over\langle
a^{2}\rangle}\right)\left({m_{\rm g}\over\rho_{\rm g}v_{\rm
th}}\right)\left({\rho_{\rm d}\over\rho_{\rm g}}\right)^{-1},$ (7)
$\displaystyle\simeq$ $\displaystyle 30.4~{}{\rm yr}\left(\rho_{s}\over
3~{}{\rm g~{}cm^{-3}}\right)\left({\langle a^{3}\rangle/\langle
a^{2}\rangle\over 0.1~{}\mu{\rm m}}\right)\left({T_{\rm g}\over 20~{}{\rm
K}}\right)^{-1/2}$ $\displaystyle\qquad\times\left({\rho_{\rm g}\over
10^{-15}~{}{\rm g~{}cm}^{-3}}\right)^{-1}\left({\rho_{\rm d}/\rho_{\rm g}\over
0.01}\right)^{-1},$

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_{\rm th}$ is the thermal velocity of the gas, and $\rho_{d}\,/\,\rho_{g}$
is the local dust-to-gas ratio. We note that $t_{\rm coll}$ therefore scales
with the local mean grain size, the local gas-to-dust ratio and inversely with
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_{\rm s}=3~{}{\rm g~{}cm^{-3}}$, $\langle
a^{3}\rangle\,/\,\langle a^{2}\rangle=0.1\micron$, $\rho_{\rm
g}=10^{-15}~{}{\rm g~{}cm^{-3}}$ (an estimate for one scale height above the
midplane assuming a surface density $\Sigma_{\rm g}=1~{}{\rm g~{}cm^{-2}}$ and
a gas scale height of 20 au, broadly consistent with the model of Calahan et
al., 2021), $T_{\rm g}=20$ K, and $\rho_{d}\,/\,\rho_{g}=0.01$, the
collisional timescale is $t_{\rm coll}\simeq 30.4~{}{\rm years}$. This
representative value is several times larger than the $t_{\rm coll}$ inferred
from the observed shadow offset in TW Hya (bottom panels of Fig 5 and 5). Such
a discrepancy could easily be accounted for with variations of a factor of a
few for all the properties described in Eqn. 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.,
Fig. 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\farcs 4$),
before the CO emission was no longer detectable. The difference between these
archival observations and the ones presented here are that the previous
observations were at slightly higher frequencies, observing the $J=3-2$
transition with $E_{u}=33.2~{}{\rm K}$, as opposed to the lower energy
$E_{u}=16.6~{}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 a lower energy transitions will probe deeper layers in 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_{\rm coll}$. When
$t_{\rm coll}$ increases to become a substantial fraction of $t_{\rm 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_{\rm 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_{\rm coll}(r)<{\Delta\phi_{s,\,{\rm dust}}(r)\over|\Omega_{s}-\Omega_{\rm
kep}(r)|},$ (8)

where $\Delta\phi_{s,\,{\rm 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\fdg 9~{}{\rm yr^{-1}}$ (as any smaller $\Omega_{s}$ will just
increase the limits of $t_{\rm coll}$). As $|\Omega_{s}|\gg|\Omega_{\rm kep}|$
at these radii, it is irrelevant if the disk is rotating clockwise or
anticlockwise, as the relative motion is dominated by the movement of the
shadow. We find $t_{\rm coll}\lesssim 8~{}{\rm years}$ is required for a gas
to feel the effects of the shadow. A shadow crossing time of ${\sim}~{}8$
years sits squarely in the range of $t_{\rm 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_{\rm
coll}$ is close to this limit, and $J=3-2$ probes an elevated layer where
$t_{\rm coll}>8$ years 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\arcsec$. As the
surface density drops towards 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_{\rm coll}$ values until a threshold is reached (Eqn. 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 exists.

##  7 Summary

In this paper we presented 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 reanalyzed. It
was found 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° – 60° 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_{\rm coll}$ of between 1 and 10 years.

Although the HST data is 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_{\rm coll}$, ranging between 5 and
20 years.

The inferred $t_{\rm 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_{\rm coll}$ is
larger due to the lower dust densities.

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 Union’s Horizon 2020
research and innovation programme (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.

## References

  * Akiyama et al. (2015) Akiyama, E., Muto, T., Kusakabe, N., et al. 2015, ApJ, 802, L17, doi: 10.1088/2041-8205/802/2/L17
  * Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40, doi: 10.3847/2041-8205/820/2/L40
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Apai et al. (2004) Apai, D., Pascucci, I., Brandner, W., et al. 2004, A&A, 415, 671, doi: 10.1051/0004-6361:20034549
  * Bae et al. (2021) Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56, doi: 10.3847/1538-4357/abe45e
  * Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644, doi: 10.1038/nature11805
  * Burke & Hollenbach (1983) Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223, doi: 10.1086/160667
  * Calahan et al. (2021) Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8, doi: 10.3847/1538-4357/abd255
  * Canta et al. (2021) Canta, A., Teague, R., Le Gal, R., & Öberg, K. I. 2021, ApJ, 922, 62, doi: 10.3847/1538-4357/ac23da
  * Casassus et al. (2019) Casassus, S., Pérez, S., Osses, A., & Marino, S. 2019, MNRAS, 486, L58, doi: 10.1093/mnrasl/slz059
  * Cleeves et al. (2015) Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204
  * Czekala et al. (2021) Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2, doi: 10.3847/1538-4365/ac1430
  * Dartois et al. (2003) Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773, doi: 10.1051/0004-6361:20021638
  * Debes et al. (2016) Debes, J. H., Jang-Condell, H., & Schneider, G. 2016, ApJ, 819, L1, doi: 10.3847/2041-8205/819/1/L1
  * Debes et al. (2013) Debes, J. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45, doi: 10.1088/0004-637X/771/1/45
  * Debes et al. (2017) Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205, doi: 10.3847/1538-4357/835/2/205
  * Dullemond et al. (2020) Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137, doi: 10.1051/0004-6361/201936438
  * Facchini et al. (2017) Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16, doi: 10.1051/0004-6361/201630329
  * Facchini et al. (2018) Facchini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459, doi: 10.1093/mnras/stx2523
  * Flaherty et al. (2018) Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117, doi: 10.3847/1538-4357/aab615
  * Foreman-Mackey et al. (2019) Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, The Journal of Open Source Software, 4, 1864, doi: 10.21105/joss.01864
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gorti et al. (2011) Gorti, U., Hollenbach, D., Najita, J., & Pascucci, I. 2011, ApJ, 735, 90, doi: 10.1088/0004-637X/735/2/90
  * Huang et al. (2018) Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018, ApJ, 852, 122, doi: 10.3847/1538-4357/aaa1e7
  * Hughes et al. (2011) Hughes, A. M., Wilner, D. J., Andrews, S. M., Qi, C., & Hogerheijde, M. R. 2011, ApJ, 727, 85, doi: 10.1088/0004-637X/727/2/85
  * Jorsater & van Moorsel (1995) Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037, doi: 10.1086/117668
  * Kama et al. (2016) Kama, M., Pinilla, P., & Heays, A. N. 2016, A&A, 593, L20, doi: 10.1051/0004-6361/201628924
  * Krist et al. (2000) Krist, J. E., Stapelfeldt, K. R., Ménard, F., Padgett, D. L., & Burrows, C. J. 2000, ApJ, 538, 793, doi: 10.1086/309170
  * Law et al. (2021) Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4, doi: 10.3847/1538-4365/ac1439
  * Loomis et al. (2018) Loomis, R. A., Cleeves, L. I., Öberg, K. I., et al. 2018, ApJ, 859, 131, doi: 10.3847/1538-4357/aac169
  * Macías et al. (2021) Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33, doi: 10.1051/0004-6361/202039812
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Muley et al. (2021) Muley, D., Dong, R., & Fung, J. 2021, AJ, 162, 129, doi: 10.3847/1538-3881/ac141f
  * Nealon et al. (2019) Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951, doi: 10.1093/mnras/stz346
  * Nealon et al. (2020) Nealon, R., Price, D. J., & Pinte, C. 2020, MNRAS, 493, L143, doi: 10.1093/mnrasl/slaa026
  * Qi et al. (2013) Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Science, 341, 630, doi: 10.1126/science.1239560
  * Rapson et al. (2015) Rapson, V. A., Kastner, J. H., Millar-Blanchaer, M. A., & Dong, R. 2015, ApJ, 815, L26, doi: 10.1088/2041-8205/815/2/L26
  * Roberge et al. (2005) Roberge, A., Weinberger, A. J., & Malumuth, E. M. 2005, ApJ, 622, 1171, doi: 10.1086/427974
  * Schwarz et al. (2016) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91, doi: 10.3847/0004-637X/823/2/91
  * Teague (2019a) Teague, R. 2019a, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632
  * Teague (2019b) —. 2019b, The Journal of Open Source Software, 4, 1220, doi: 10.21105/joss.01220
  * Teague (2020) Teague, R. 2020, richteague/keplerian_mask: Initial Release, 1.0, Zenodo, Zenodo, doi: 10.5281/zenodo.4321137
  * Teague et al. (2019a) Teague, R., Bae, J., & Bergin, E. A. 2019a, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0
  * Teague et al. (2019b) Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019b, ApJ, 884, L56, doi: 10.3847/2041-8213/ab4a83
  * Teague & Foreman-Mackey (2018) Teague, R., & Foreman-Mackey, D. 2018, Research Notes of the American Astronomical Society, 2, 173, doi: 10.3847/2515-5172/aae265
  * Teague et al. (2016) Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49, doi: 10.1051/0004-6361/201628550
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132, doi: 10.3847/1538-4357/aa5d68
  * Walsh et al. (2016) Walsh, C., Loomis, R. A., Öberg, K. I., et al. 2016, ApJ, 823, L10, doi: 10.3847/2041-8205/823/1/L10
  * Weinberger et al. (2002) Weinberger, A. J., Becklin, E. E., Schneider, G., et al. 2002, ApJ, 566, 409, doi: 10.1086/338076
  * Young et al. (2004) Young, K. E., Lee, J.-E., Evans, Neal J., I., Goldsmith, P. F., & Doty, S. D. 2004, ApJ, 614, 252, doi: 10.1086/423609
  * Zhang et al. (2017) Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nature Astronomy, 1, 0130, doi: 10.1038/s41550-017-0130