Gas and Dust Shadows in the TW Hydrae Disk

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 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 ~100 deg. 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}$ ~ 10 years, depending on whether a straight or a 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}$.


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. • 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 Weinberger et al. 2002;Apai et al. 2004;Roberge et al. 2005;Debes et al. 2013Debes et al. , 2016Akiyama 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 ) 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 50 au (Roberge et al. 2005), appears to move across the face of the disk at a rate of 22. • 7 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. (2019Nealon et al. ( , 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 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.

OBSERVATIONS
Observations were acquired as part of project 2018.A.00021.S (PI: Teague), designed for high spatial and spectral resolution observation of 12 CO J = 2 − 1, CN N = 2 − 1 and CS J = 5 − 4. The correlator was set up to cover the 12 CO and CS line with a channel spacing of 15.3 kHz (resulting in a 30.5 kHz resolution, equivalent to 40 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 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 CH 2 CN 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.

Calibration
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 selfcalibrated following the procedure used in the DSHARP program . 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.

Imaging
Briggs weighting and a channel spacing of 40 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. 19 × 0. 17 (89 • ). A Keplerian mask was generated, 1 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 −1 . Integrating over the Keplerian mask we recover an integrated intensities of 18.31 ± 0.03 Jy km s −1 . This uncertainty does not include a systematic uncertainty of ∼ 10% associated with the flux calibration of the data.
The package bettermoments 2 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 Fig. 1b. The analysis of additional moment maps is the focus of a companion paper (Teague et al., in prep.).

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: ∆RA = 20.0 ± 0.1 mas and ∆Dec = 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 was 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 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 largescale 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 misspecified center, the same residuals were calculated for a range of offsets in both the x-and y-direction of ±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, where I(r) is the azimuthally averaged integrated intensity, δI(r) is the radially dependent amplitude of the azimuthal variation in integrated intensity and φ s, gas (r) is the radially dependent position angle of the shadow (note that unlike Debes et al. 2017 we have included a π phase offset such that φ 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 φ s, 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 {δI, φ s, gas } for each annulus, using 1024 walkers taking 20,000 steps, with the first 10,000 being discarded for burnin. Uniform priors were adopted for both parameters, with δI assumed to be positive and less than 1 mJy beam −1 km s −1 , while 0 • ≤ φ s, gas < 360 • , 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. 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 ≈ 0.4 mJy beam −1 km s −1 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 Fig. 1a). 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.

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. • 7 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 • (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.

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 where Ω s is the orbital frequency of the shadow, c the speed of light, and (r rim , φ rim (t)) describing the cylindri-cal coordinates of the shadowing rim. Adopting Ω s = 22. • 7 yr −1 , the apparent orbital frequency of the shadow, we find dφ s, dust /dr ≈ 4 × 10 −4 deg au −1 , 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 * /r 3 rim , and solved for r rim . Determining r rim for each of the six spatially resolved epochs yielded values of r rim ranging between r rim = 0.04 ± 0.01 au for 2004 to r rim = 0.8 ± 0.6 au for 2015. Such a vast variation in r rim is hard to reconcile with a model of a precessing inner disk (e.g. Facchini et al. 2018;Nealon et al. 2019Nealon et al. , 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.

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. 3 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, whereΩ s is the angular frequency of the shadow and δt is a phase offset. As in Section 3, emcee was used to sample the posterior distributions of {Ω s , δt} using 1024 walkers taking 20,000 steps, with the first 10,000 being discarded for burnin. 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Ω s = 22. • 9 ± 0. • 1 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. 2019Nealon et al. , 2020. To check the suitability of this model, we fit the data with the form, where φ s, dust is the average position angle of the shadow, δφ s, dust is the amplitude of the oscillation, τ is the period of the oscillation and δ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 δφ s, dust = 55 • ± 2 • , φ s, dust = 239 • ± 2 • and τ = 15.8 ± 0.1 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 δφ s, dust = 57 • ± 2 • , φ s, dust = 232 • ± 3 • and τ = 15.8 ± 0.1 yr. The similarity in inferred parameters suggests that the original datasets 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 φ s, dust = 288 • ± 37 • (where the quoted uncertainty describes the uncertainty-weighted standard deviation of sample), while the predicted locations werē φ s, dust = 319 • ± 1 • andφ s, dust = 288 • ± 2 • , 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 φ s, dust = 298 • ± 20 • , for which the orbital model predicts a location ofφ s, dust = 350 • ± 2 • while the oscillatory model predictsφ s, dust = 281 • ±3 • , 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 unam-biguously distinguish between these scenarios and, as such, both models are used in the remainder of this work.

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 where c V is the specific heat capacity of the gas, k B is the Boltzmann constant,α T ≈ 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 ≈ 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 the 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 shadows of ∼ 60 • − 0 • between 3. 1 and 3. 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  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 desaturated version of Fig. 1c. Note that disk is rotating in a clockwise direction.
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 where Ω s − Ω 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 Ω kep varies radially, Ω 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 φ), in the opposite direction to the rotation of the disk such that ∂φ s, dust /∂t is positive and Ω kep is negative.
As Ω s is considerably larger than Ω 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. Adopting the same stellar mass as used for the Keplerian CLEAN mask in Section 2 and the instantaneous an- , 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 Eqn. 2 (i.e., an arithmetic spiral), and linearly interpolated between the best-fit models of the 2016 and 2021 epochs for φ s, dust to model the curvature of the shadow during the ALMA observations. This results in a morphology with a radial dependence of dφ s, dust /dr ≈ 0. • 29 au −1 , such that average offset between shadows are increased by approximately 70 • , as shown in Fig. 5c. 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 Fig. 5d). 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.

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 3 ), likely related to a shadow observed in scattered light at smaller radii (r 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.

Collisional Timescale
Following Bae et al. (2021), the collisional timescale for gas molecules is given by 30.4 yr ρ s 3 g cm −3 a 3 / a 2 0.1 µm where ρ s is the bulk density of the dust grains, a 3 / a 2 is the mean dust grain size, m g is the mean mass of gas molecules, ρ g is the local gas volume density, v th is the thermal velocity of the gas, and ρ d / ρ g is the local dust-togas ratio. We note that t 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., ρ s = 3 g cm −3 , a 3 / a 2 = 0.1µm, ρ g = 10 −15 g cm −3 (an estimate for one scale height above the midplane assuming a surface density Σ g = 1 g cm −2 and a gas scale height of 20 au, broadly consistent with the model of Calahan et al. 2021), T g = 20 K, and ρ d / ρ g = 0.01, the collisional timescale is t coll 30.4 years. This representative value is several times larger than the t 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, ρ d /ρ 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).

Excitation Conditions
The shadows described in this work were not observed in previously published high resolution ALMA data of CO emission 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 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 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 coll . When t coll increases to become a substantial fraction of t 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 coll than the shadow-crossing time this could explain why shadows were not seen, but spirals were (which Bae et al. 2021 andMuley 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 where ∆φ s, 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 ∆φ s 180 • for the shadow and the maximum angular velocity of the shadow, Ω s = 22. • 9 yr −1 (as any smaller Ω s will just increase the limits of t coll ). As |Ω s | |Ω 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 coll 8 years is required for a gas to feel the effects of the shadow. A shadow crossing time of ∼ 8 years sits squarely in the range of t 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 coll is close to this limit, and J = 3 − 2 probes an elevated layer where t 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 ∼ 3 . 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 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.

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 ∼ 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 nonnegligible collisional timescale between the gas molecules and dust grains. Given the uncertainty in the angular velocity of the shadow, this corresponds to a t 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 coll , ranging between 5 and 20 years.
The inferred t 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 coll is larger due to the lower dust densities.