The Surprising Evolution of the Shadow on the TW Hya Disk^*

TW Hya is a nearby, well-studied face-on protoplanetary disk (d = 60 pc, i ∼ 5°; Qi et al. 2004; Huang et al. 2018; Gaia Collaboration et al. 2021). It has been observed from X-ray to radio wavelengths. Estimates of TW Hya's stellar mass generally lie in the range 0.6–0.8 M_⊙, with recent (dynamical) estimates favoring the upper end of this range (e.g., Teague et al. 2022). Given its age (∼10 Myr), this makes the TW Hya disk an interesting test case for planet formation in a large, relatively isolated, gas-rich disk orbiting a star somewhat less massive than the young Sun but more massive than transiting young M-dwarf systems such as AU Mic (Sokal et al. 2018; Teague et al. 2019). The Hubble Space Telescope (HST) has imaged the disk numerous times, including with WFPC2 (inner working angle, IWA ≈ 07; Krist et al. 2000), NICMOS near-IR (NIR) F110W and F160W (IWA ≈ 04; Weinberger et al. 2002), and Space Telescope Imaging Spectrograph (STIS) visible-light imaging/spectroscopy (IWA ≈ 02; Roberge et al. 2005; Debes et al. 2017, hereafter D17). The disk is polarized, allowing visible and NIR images (Akiyama et al. 2015; Boekel et al. 2017, hereafter VB17; Poteet et al. 2018). In the submillimeter, dust continuum emission extends from 3 to ∼100 au (Andrews et al. 2016; Huang et al. 2018; Ilee et al. 2022), with multiple gaps or rings present that are close in radius to structures seen in scattered light and consistent with active planet formation (Jang-Condell & Turner 2012; Andrews et al. 2016; Ueda et al. 2020; Macías et al. 2021, and references therein). CO gas emission extends to ∼200 au, and scattered light from small dust extends as far out as 452 au (Huang et al. 2018; D17).

The TW Hya disk has an azimuthal surface brightness (ASB) asymmetry between 30 and 150 au (Roberge et al. 2005; Debes et al. 2013). The asymmetry has been attributed to inclination effects for flared disks (Debes et al. 2013), the forward-scattering properties of dust (Roberge et al. 2005), or the presence of a warp (Roberge et al. 2005; Rosenfeld et al. 2012). Images of the asymmetry with STIS and NICMOS over six epochs between 0.6 and 2 μm show that the asymmetry is not stationary; the apparent changes in disk illumination can be modeled with a rotating shadow with a constant angular velocity of 22.7° yr⁻¹ (P = 15.9 yr). Assuming the shadow arises from an obstructing feature moving with Keplerian velocity constrains the location of the originating structure to $\lt 5.6(\tfrac{{M}_{\star }}{0.7{M}_{\odot }})$ au (D17). Atacama Large Millimeter/submillimeter Array (ALMA) CO maps of the outer disk also show the shadow (Teague et al. 2022). The position of the shadow in CO taken in 2019 is smaller than expected based on the predictions of D17.

The TW Hya shadow is one of a handful of other disks that show evidence of inner occulting structures, often interpreted as misaligned inner disks casting narrow or broad shadows depending on the inclination difference with the outer disk (e.g., PDS 66, HD 34700, HD 16192, HD 142527, HD 139614, and RXJ1604.3-2130A; Marino et al. 2015; Wolff et al. 2016; Bertrang et al. 2018; Pinilla et al. 2018; Muro-Arena et al. 2020; Uyama et al. 2020). Some of these disks also show evidence for time variability, but robustly sampling the variability timescales is difficult for high-contrast imaging targets. The robustness of the shadow detection and coverage of its characteristic timescale of motion make the shadow on TW Hya's disk one of the best-characterized shadows to date that may be caused by active planet formation.

Nealon et al. (2019) investigated the behavior of shadows due to misaligned disks that precess and warps associated with planetary companions orbiting the host star, in order to understand TW Hya's shadow. They found that massive planets with inclinations relative to the outer disk of a few degrees or more were sufficient to create shadows of the right amplitude in the outer disk either through a warp exterior to the planet's orbit or from a misaligned disk interior to the planet. This work motivated a multicycle HST program with the STIS coronagraph in order to disentangle these two origins. TW Hya has only been observed for about half the observed periodicity, and a misaligned, precessing disk will continue to show asymmetries in the disk for a whole precession period, while features associated with a misaligned planet will disappear for half an orbital period. It is possible that both structures are present and contributing to the shadow behavior.

We report here the first-epoch results from our monitoring program of TW Hya's disk to further track the evolution of the asymmetry and putative shadow. In Section 2 we describe the new STIS observations. In Section 3 we report our findings of the latest epoch, which show a marked departure in the behavior of the shadow. In Section 4 we discuss the possible explanation of this evolution. Finally, we present our conclusions in Section 5.

TW Hya was observed on 2021 June 7 with the STIS instrument (Woodgate et al. 1998), using the charge-coupled device (CCD) and 50CORON aperture as part of GO#16228 (PI: Debes). The same point-spread function (PSF) reference for TW Hya, HD 85512, was used as in Roberge et al. (2005) and D17. TW Hya was observed for two consecutive orbits with spacecraft orientations separated by 17°, followed by an orbit on HD 85512, and a final orbit on TW Hya separated by 32° from the first TW Hya orbit. Each star was first placed behind the BAR5 aperture location, which has a half-width of 015 with inner working angles of ∼02 (Schneider et al. 2017; Debes et al. 2019). For TW Hya, 5 × 110 s exposures were obtained per orbit behind BAR5, with sufficiently short exposure times chosen to avoid saturating the detector at the mask edge. The PSF reference was dithered perpendicular to BAR5 by ±12.5 mas to account for the 0.25 pixel target-acquisition nonrepeatability (Schneider et al. 2017). We obtained 5 × 11.2 s exposures at each dither location. Next, each star was placed behind the WEDGEA1.0 occulter to ensure sufficient signal-to-noise ratio (S/N) in the outer disk without saturating the detector. The WEDGEA1.0 exposures were 3 × 425 s per orbit with TW Hya and 22 × 32 s exposures for HD 85512. The total integration time for TW Hya was 1100 s behind BAR5 and 2550 s behind WEDGEA1.0. This procedure replicates the other STIS images obtained in 2015 and 2016 (D17).

Unfortunately, a delay in guide-star acquisition rendered the final orbit of TW Hya unusuable, since the star was not placed behind the proper apertures. The purpose of the multiple orbits of TW Hya are primarily to ensure nearly 360° coverage on the disk exterior to 05. The result of the failed orbit is a modest decrease in disk S/N and a smaller angular coverage at small angular radii, as seen in Figure 1. However, as we show below, this did not significantly impact our results when comparing to previous epochs at radii >40 au (067).

Zoom In Zoom Out Reset image size

We followed the PSF subtraction laid out in D17, which we briefly review. First, the stellar center was determined for both target and PSF reference. We then iteratively registered and scaled the PSF reference, minimizing subtraction residuals within a mask that included the diffraction spikes of the star and subtracted off the PSF reference. We note that the ratio of TW Hya's flux to HD 85512 across the STIS CCD bandpass was higher in 2021 between 0.2 and 1.0 μm by 24%–31%. Compared to 2015 and 2016 where the best-fit ratio of TW Hya's flux to HD 85512 was 0.052 and 0.055, we find that the best scaling between the two stars is now 0.0683.

TW Hya has been known to be variable by ∼20%–30% in the visible on periods close to ∼3 days (Mekkaden 1998). To investigate whether this magnitude of brightening might be expected in the STIS bandpass, we downloaded existing contemporaneous STIS/G430L and STIS/G750L spectra taken with the 52 × 2 slit of TW Hya available in MAST and estimated the count rate on the STIS detector for the 50CORON mode for each spectrum using the STIS Exposure Time Calculator (ETC). The spectra came from four epochs (2010 January 28, 2010 February 4, 2010 May 28, and 2015 April 18) as part of programs 11608 (PI: Calvet) and 13775 (PI: Espaillat). The varying flux in the blue continuum and in emission-line intensity of the different spectra account for peak-to-peak changes in predicted count rates on the STIS detector across the four epochs of 22%, in line with the variability observed. Additionally, we inspected the acquisition images for TW Hya (which used the F28X50LP filter) between 2016 and 2021. Even though the filter throughput for the acquisition images is for flux beyond 5000 Å, aperture photometry of the acquisition images showed a 22% increase in TW Hya's brightness between the two epochs. The change in the spectral energy distribution (SED) for TW Hya adds modest variability in how well the target PSF reference matches in color to TW Hya, which could result in more or less PSF subtraction residuals.

With the scalings, we recover the disk again, as seen in Figure 1. We show reductions of each individual visit to demonstrate the impact of PSF subtraction residuals. Residual features are present only in the detector frame of the observations and will appear to rotate when in a sky-oriented frame. While some modest features appear to be due to PSF residuals, they are at the level of our estimated noise per pixel and do not impact our measurements significantly for the analysis described in Section 3.

As noted in D17 and VB17, when one removes the 1/R² dependence of disk illumination from the star, a gap is present at ∼80 au, with an inner ring peaking at ∼30 au, with another gap interior to that, as shown in the 2016 image in Figure 2. There is an ASB asymmetry, which we interpret as due to a shadow (D17). In past epochs, this shadow asymmetry was well described by a smooth sinusoidal curve. The position angles (PAs) of the peak and trough of the asymmetry have, to date, moved with constant angular velocity across several wavelengths in the visible and NIR. Using the 15.9 yr predicted period from D17 and setting t_o to equal the first observation of the TW Hya disk, we find that roughly half of the predicted orbit for the shadow has not been covered, and the 2021 June observations mark a similar phase of the orbit to archival NICMOS polarimetric observations (Poteet et al. 2018) and archival medium-band NICMOS total-intensity observations of the disk (D17) taken in 2004 and 2005, both of which are of lower S/N compared to the observations with STIS.

We created an azimuthal average surface brightness profile of the 2021 data set, converted it into physical flux units per pixel, and compared the resulting fluxes with azimuthal averages of the 2016 observations of the TW Hya disk (see Figure 4). Due to the broad bandpass of the 50CORON mode, we estimate the conversion from counts per second on the detector to surface brightness by using the stsynphot package and the average spectrum of TW Hya obtained with STIS. We find that the conversion is 1 mJy arcsec⁻² = 3.54 DN s⁻¹. We find that the flux profile of TW Hya is slightly different in 2021 compared to 2016, primarily fainter by 20% interior to 1'' as well as beyond 15.

Unlike the observations taken in 2016, the inner ring shows two narrow dark lanes at PAs of ∼50° and ∼180°. The peak of the azimuthal asymmetry has moved counterclockwise toward due north as predicted by the shadow rotation model presented in D17. To quantify these changes, we replicate the methods of D17 by measuring the ASB of the disk as a function of PA and radius for the 2021 epoch and compare those to the previously published profiles in 2016.

To obtain azimuthal asymmetries, we divided each ASB profile for a given radius by the mean surface brightness at that radius to allow for unbiased radial averaging. We chose angular radii that matched the measurements made in D17, and converted to physical differences with the new Gaia Early Data Release 3 (EDR3) parallax (Gaia Collaboration et al. 2016, 2021): seven radial locations centered at 046 (28 au), 066 (40 au), 089 (54 au), 114 (68 au), 147 (88 au), 183 (110 au), and 236 (142 au), with widths of 02 (12 au), 02 (12 au), 025 (15 au), 025 (15 au), 041 (25 au), 030 (18 au), and 056 (34 au), respectively. The locations are within and near the two gaps seen in scattered light at 22 and 88 au with sufficient S/N. We estimated the rms uncertainty per azimuthal bin by calculating both the average uncertainty per pixel as estimated by the total counts within each pixel or by the standard deviation of counts within an azimuthal bin, whichever was larger. At these radii, the profiles in 2016 beyond 40 au have a median S/N of 19 per azimuthal bin, while for 2021 the median S/N per azimuthal bin is 18.

Zoom In Zoom Out Reset image size

Figure 3 shows the results for the 2016 and 2021 epochs. The simple cosine fit used in D17 is no longer sufficient to explain the observed ASB profile in 2021, thanks to the apparent presence of two shadow features. Figure 3 shows that the two shadows diminish in depth at larger radii, such that they are no longer distinguishable as two components and that the peak surface brightness increases in PA with radius. We thus fit a cosine at 141 au and the PA of the peak surface brightness is 330° ± 10°, in tension (∼2.3σ) with the 2005 measurement from the NICMOS medium-band data at a similar phase of the predicted motion reported in D17. We note that the 2004 NICMOS total-intensity data derived from the polarizing filters in D17 and Poteet et al. (2018; P19) showed a similarly lower-than-expected PA for the given period.

Zoom In Zoom Out Reset image size

There are three possibilities that can explain the new observations. First, there are new features caused by the recent appearance of additional shadowing structures unrelated to what has been seen before, such as new dust clumps. Second, warps due to companions that cast shadows can appear when the planet is on the near side of the disk from the perspective of the observer (Nealon et al. 2019). This is possible since we are seeing the phenomena at a PA that is close to the major axis of the disk at 152° as measured by CO data (Huang et al. 2018). Third, it is possible that these two shadows arise from two mutually inclined disks that have precession timescales which are slightly different and casting overlapping shadows, thus mimicking a single shadow over the previously observed epochs.

The ASB profiles can be fit by combinations of shadows and cosine functions. At present it is difficult to find a unique solution between these possibilities. We therefore focus on the simplest case: that there are two shadows which previously overlapped but are moving at different angular velocities and partially separated in 2021. Thus, we can approximate the shape of the shadows as a linear combination of Gaussian curves.

In order to fit the two shadow features simultaneously, we use CURVEFIT.pro, an IDL nonlinear least-squares fitting routine assuming two Gaussian shadows, with constant offset (C), shadow depth (D_i), shadow PA (θ_o,i), and shadow angular width (σ_i):

$$\begin{eqnarray}&&F(\theta )=C+\displaystyle \sum _{i=1}^{2}{D}_{{\rm{i}}}\exp \displaystyle \frac{-{\left(\theta -{\theta }_{{\rm{o}},{\rm{i}}}\right)}^{2}}{2{\sigma }_{{\rm{i}}}^{2}},\end{eqnarray} \tag{ 1 }$$

where we convert the shadow width to an angular FWHM, FWHM_i ≈ 2.3548σ_i. The fits are shown in Figure 3 and show no significant residuals.

Zoom In Zoom Out Reset image size

Figure 5 shows the PAs, widths, and depths of the shadows as a function of radius. Shadow A has a median PA of 49°, while Shadow B has a median PA of 175°. The shadows show a slight increase in PA in the counterclockwise direction as radius increases, which may be indicative of warping in the disk's surface.

Zoom In Zoom Out Reset image size

The FWHMs of the shadows are broad: the median FWHM of shadow A is 87°, shadow B 93°. These widths can explain why the shadow might have appeared to be singular previously if the separation of the two shadows was less than their FWHMs. Finally, the depth of the shadow changes with radius, starting at depths of near 80% and decreasing to 20% by 140 au.

We review the following characteristics and evolution of the disk shadow in order to provide a summary of the current observational constraints. We will restrict ourselves to those behaviors that are seen in the STIS data, since this provides a set of constraints that cannot be explained by changes in behavior of the shadow with wavelength, instrumental setup/systematics, or disk polarization.

From 2000 to 2016 the shadow covered 180° of azimuth and the motion of both edges of the shadow were consistent with constant angular motion. The shadow depth was constant over this time period and with radius. The PA of the shadow at r < 50 au did not always match that of outer parts of the disk. The disk shadow detected prior to 2016 may have evolved into two shadows by 2021. The behavior of the shadow between 2000, 2015, and 2016 suggested that the shadow rotated counterclockwise with a period of ∼15.9 yr. The behavior of the shadow in 2021 may mean that either rotation or oscillatory behavior of the shadow is possible, but the overall period still seems to be ∼15 yr, limited somewhat by sparse time sampling of the shadow's motion with STIS-only epochs. In 2021, there is a trend where the depth of the shadows changes quasi-monotonically with radius, with the shadow becoming less pronounced further out in the disk. This behavior is not as pronounced in other STIS epochs.

In order to assign a quasi-physical explanation to this new shadow behavior, we begin in Section 4.1 to model the disk without any shadowing structures as a baseline disk model to compare to two different scenarios: a single inclined ring casting a shadow, and two mutually inclined rings casting shadows. In Section 4.2 we attempt to fit the data with one inclined ring. This is because a tilted ring is physically motivated (e.g., D17), conserves angular momentum, in some cases has been observed (Marino et al. 2015), and has been inferred in other disks (Muro-Arena et al. 2020). We then use the rapid radiative-transfer code MCFOST (Pinte et al. 2006, 2009) to model the shadow cast by an inclined ring to try and change the location, shape and inclination of this ring to match the flux profiles in Figure 3. Ultimately, we show the best-fitting inclined inner ring cannot fit the observed ASB and investigate the impact of adding a second inclined ring in Section 4.3. There we find a good match to both the SED of the disk and the visible scattered-light images, which suggests that this scenario is the most likely explanation for the behavior we observe.

4.1. A Fiducial Model for the TW Hya Disk

The first step is to create a disk model that is broadly consistent with both the visible total-intensity scattered-light morphology and the SED measured for TW Hya without any shadows. To that end, we recreated the model proposed by VB17 with MCFOST. VB17 was focused on simultaneously fitting the polarized-intensity radial surface brightness profile of the disk NIR wavelengths with the TW Hya SED. This model was also derived from work done by Menu et al. (2014), which was a model focused on self-consistently fitting multiwavelength interferometry data with the TW Hya SED. Both models relied on a one-dimensional axisymmetric disk structure with a vertical scale height dictated by hydrostatic equilibrium. Both models assumed that the inner gas disk was truncated with a rounded inner taper and that the largest dust grains were decoupled from the gas radial density profile. VB17 additionally added an outer gas density taper at a larger radius than Menu et al. (2014), used the ALMA 870 μm radial surface brightness profile (Andrews et al. 2016) as a proxy for the large dust radial profile, and added gaps in the form of density decrements at particular radii to reproduce the polarized-intensity radial scattered-light surface brightness profiles observed with SPHERE. For the purposes of our models, we assumed that the outer disk has an inclination of 0° and a PA on the sky of 155°. While the true inclination of the outer disk is likely closer to 5°, the impact to both the scattered-light images and SED are negligible.

While we consider a very similar disk structure, we make some alternate choices that appear to retain the basic observational features of the TW Hya disk while making the model more amenable to calculating the radiative transfer of tilted inner rings. We are not trying to make a disk model that matches all available observations but a plausible model that can be used to investigate the impact of tilted inner rings on the observed surface brightness. We leave a deeper investigation into the degeneracies and constraints imposed by disk shadows for later work. In general, we find decent agreement with those of VB17, though we point out differences where they occur.

4.1.1. Radial Surface Density of Gas and Small Dust

The previous two models of the TW Hya disk calculated an analytical radial surface density profile for the gas while assuming a gas-to-dust ratio for the small dust grains and calculated a self-consistent scale height for the disk as a function of radius. Our approach is to calculate a fully analytical gas density distribution with the following form:

$$\begin{eqnarray}&&\rho (r,z)={{\rm{\Gamma }}}_{\mathrm{in}}(r)\displaystyle \frac{{\rm{\Sigma }}(r)}{\sqrt{2\pi }H(r)}f(r){{\rm{\Gamma }}}_{{\text{}}{out}}{\int }_{-\infty }^{\infty }\exp \left(-\displaystyle \frac{{z}^{2}}{2H{\left(r\right)}^{2}}\right){dz},\end{eqnarray} \tag{ 2 }$$

where Γ_in(r) is the interior density taper, Σ(r) is the surface density profile, f(r) is a multiplicative combination of density gaps to mimic the gaps observed in various scattered-light data sets, Γ_out is the exterior density taper, and H(r) is the vertical scale height of the gas (VB17, and references therein).

The interior taper replicates the format in VB17, and is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{in}}(r\lt {r}_{\mathrm{in}})=\exp \left[{\left(\displaystyle \frac{1-r/{r}_{\mathrm{in}}}{w}\right)}^{3}\right].\end{eqnarray} \tag{ 3 }$$

The surface density profile is a power law:

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={\left(\displaystyle \frac{r}{{r}_{\mathrm{in}}}\right)}^{-p},\end{eqnarray} \tag{ 4 }$$

where we choose p = −0.7, and r_in is the location where the interior taper takes effect.

The scattered-light profile of TW Hya clearly shows the evidence of gaps or other perturbations in the upper layers of the disk, some of which correspond to features also seen in the submillimeter where the disk is optically thin beyond a few 10 s of au (Debes et al. 2013; Andrews et al. 2016; VB17). We reproduce these features by including Gaussian density perturbations following VB17:

$$\begin{eqnarray}&&{f}_{i}=1-{d}_{i}\exp \left[-\displaystyle \frac{{\left(r-{r}_{c,i}\right)}^{2}}{2{\sigma }_{\mathrm{in}},{i}^{2}}\right]\,\mathrm{for}\ r\lt {r}_{c,i},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{f}_{i}=1-{d}_{i}\exp \left[-\displaystyle \frac{{\left(r-{r}_{c,i}\right)}^{2}}{2{\sigma }_{\mathrm{out},i}^{2}}\right]\,\mathrm{for}\ r\gt {r}_{c,i}.\end{eqnarray} \tag{ 6 }$$

While we assume the perturbations are in density, we do not include similar perturbations in scale height, which might also occur in a system undergoing gap formation (Bi et al. 2021).

We include a taper to the outer part of the gas density to match the sharp drop in surface brightness seen in scattered light and in deep observations of the outer gas disk in the submillimeter (D17; Ilee et al. 2022). We assume the following functional form for the taper that occurs at r > r_out:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{out}}=\exp \left[\displaystyle \frac{-{\left(r-{r}_{\mathrm{out}}\right)}^{2}}{2{\sigma }_{\mathrm{out}}}\right]\mathrm{for}\ r\gt {r}_{\mathrm{out}}.\end{eqnarray} \tag{ 7 }$$

Finally, we assume an analytical scale height such that $\tfrac{H}{r}$ = 0.05 at 0.5 au and that the disk has a flared geometry:

$$\begin{eqnarray}&&H(r)=0.05{\left(\displaystyle \frac{r}{0.4\mathrm{au}}\right)}^{1+q},\end{eqnarray} \tag{ 8 }$$

where q = 0.15, shallower than the flaring expected for a purely isothermal disk (Chiang & Goldreich 1997). We choose this analytical form to best match the SED and the surface brightness of the observed disk. To realize the density distribution we create a cylindrical grid in the density with 4.8 × 10⁵ cells logarithmically spaced in radius.

The detailed parameters used for our fiducial model are listed in Table 1 and 2. We compare these values to those used in VB17 and find that they are quite similar. The biggest differences occur for the location of the gap in the outer disk listed in Table 2—VB17 assume a gap at 94 au while we find a gap at 85 au better reproduces the STIS azimuthally averaged SB profile. We additionally require different widths to the gaps in Table 2 compared to those assumed in VB17.

Table 1. Table of Parameters Used to Describe the Density Distribution of Gas in the TW Hya Disk

Parameter	VB17	Fiducial	One Ring	Two Rings
p	−0.75	−0.70	−0.70	−0.70
w	0.45 au	0.50 au	0.50 au	0.50 au
rin	2.7 au	3 au	3 au	3 au
q	⋯ a	0.15	0.15	0.15
rout	104 au	100 au	100 au	100 au
σout	50 au	50 au	50 au	50 au

Note.

^a q was not a free parameter in VB17's model but was calculated via hydrostatic equilibrium.

Download table as:  ASCIITypeset image

Table 2. Table of Parameters Used to Describe the Density Gaps in the TW Hya Disk

Parameter	VB17 a	Fiducial	One Ring	Two Rings
d1	0.44	0.44	0.44	0.42
rc,1	94 au	85 au	85 au	85 au
σin, 1	16.7 au	12 au	12 au	10 au
σout,1	17.8 au	30 au	30 au	30 au
d2	0.61	0.65	0.65	0.9
rc,2	23 au	24 au	24 au	24 au
σin, 2	4.2 au	8 au	8 au	10 au
σout,2	20 au	8 au	8 au	8 au
d3	0.81	0.81	0.81	0.92
rc,3	6.7 au	6.7 au	6.7 au	6.7 au
σin, 3	3.3 au	3.3 au	3.3 au	3.3 au
σout,3	14.4 au	18.7 au	18.7 au	18 au

Note.

^a Corrected for Gaia DR3 distance.

Download table as:  ASCIITypeset image

4.1.2. Stellar Parameters

The fundamental stellar parameters of TW Hya cover a wide range in mass and T_eff, primarily because of the intrinsic difficulty of determining the parameters of young stars. For the purposes of the fiducial model, we assume similar values to those assumed in VB17, accounting for the larger distance to TW Hya determined by Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2021). In MCFOST, we chose T_eff = 4000 K and $\mathrm{log}g$ = 4.0, assuming a mass of 0.6 M_⊙. These parameters are quite close to those determined by recent studies of TW Hya's temperature and gravity (Sokal et al. 2018). They are also consistent with the inferred ${M}_{\star }\sin \,i$ from CO measurements of Keplerian velocities (Huang et al. 2018; Teague et al. 2018).

4.1.3. UV Excess

For completeness, we include a UV excess to TW Hya that fits the visible observed photometry. In reality, TW Hya's UV flux can vary by factors of 2 or more. We verified that the addition or subtraction of the UV excess in MCFOST has a minimal impact on the predicted scattered-light surface brightness and the SED.

4.1.4. Grain Properties

We initially replicated the composition of dust used in VB17, namely using a mixture of 80% amorphous magnesiumiron olivine-style silicates (MgFeSiO₄; Dorschner et al. 1995) and 20% amorphous carbon (Rouleau & Martin 1991). We note that VB17 made use of different amorphous carbon optical constants but it does not appear to greatly impact the results (Preibisch et al. 1993). We selected a Mathis et al. (1977) distribution for the dust, using a minimum grain size of 0.1 μm and maximum grain size of 10,000 μm. This effectively assumes that the large grains trace the smaller grains, which is not supported either by interferometry in the mid-IR or in the submillimeter, nor does it follow the models of VB17 or Menu et al. (2014), where large grains followed a separate density pattern. That said, since our model focused solely on reproducing a one-dimensional SED and two-dimensional scattered-light model, we assume that this choice has only a minor bearing on our results with respect to reproducing shadows. We assume that the gas-to-dust ratio in the disk is 100. The composition of dust considered by VB17 results in a surface brightness that is too faint to match the STIS observations. This can be ameliorated by increasing the scale height of the disk, but resulted in a poorer fit to the SED. After some trial and error, we elect to use a composition more similar to that of Debes et al. (2013), which is based on SED fitting by Calvet et al. (2002). For our fiducial model, we assume a mixture of 60% olivine and 40% water ice, but we note that this is likely a nonunique solution.

4.1.5. MCFOST Parameters

We input the analytical three-dimensional density profile as a grid into MCFOST and then selected 10⁷ for the number of photon packages used for the scattered-light images. Additionally, we selected a pixel size equivalent to the plate scale of a STIS pixel, we fixed the disk mass at 5 × 10⁻⁴ M_⊙, the inner radius at 0.38 au and the outer radius at 222 au to match those of VB17 accounting for the updated DR3 Gaia distance. We assumed Mie scattering for our calculations and a single wavelength at 7000 Å that mimics the expected central wavelength of the observations. The STIS CCD bandpass extends from 2000 to 10000 Å, and thus the central wavelength of an observation is dependent on the underlying SED of the source, which to zeroth order is equivalent to the SED of TW Hya, since the dust scatters light nearly constantly across visible wavelengths (Roberge et al. 2005). Additionally, we calculate a model PSF from the TinyTim software, ¹³  assuming a K7 stellar spectral type, and convolve it with the model images to compare with the STIS data.

4.1.6. Fiducial Model Fidelity

Figure 6 shows a comparison of the fiducial model SED against the observed photometry reported in Menu et al. (2014). The model SEDs for the fiducial model and the subsequent tilted-ring models we consider are virtually identical, demonstrating that such a structure has a relatively minimal impact on the one-dimensional unresolved SED. Our fiducial SED fits about as well as those of Menu et al. (2014) and VB17, though we did not attempt to optimize the fit; the reduced χ² for the fiducial fit is 3.22. Overall, the SED is close to that observed, with an overprediction of flux between 50 and 200 μm and slight underprediction of flux in the submillimeter regime. Additionally, Figure 9 shows a comparison between the observed STIS image in 2021 and the fiducial model. The agreement between both the SED and the general surface brightness of the disk is excellent. Since the fiducial disk is face-on and with no shadows, the ASB is constant.

Zoom In Zoom Out Reset image size

4.2. Shadows from Single Tilted Rings

4.3. Shadows from Two Tilted Rings

We repeat this process but consider two independently oriented tilted rings (see Figure 7). Here we choose two rings misaligned by 5.5° and 7° with respect to the outer disk plane, the first with a PA = 170° extending from 5 to 6 au and the second from 6 to 7 au with PA = 50°. These distances are consistent with the ∼16 yr timescale of the existing shadow motions assuming the motions are Keplerian and the magnitude of the inclination is chosen to roughly match the depth of the observed shadows, but other configurations are likely also possible. Investigating the degeneracies in the ring location and width is beyond the scope of this paper. If we place the rings too far out, r_c > 15 au, we cannot easily reproduce the surface brightness beyond 30 au.

Zoom In Zoom Out Reset image size

We find that the inclusion of a second ring forces slight alterations to the fiducial inner gap depths to retain a good fit to the azimuthally averaged surface brightness profile. In particular, the 24 au gap needs d = 0.9 as opposed to d = 0.65 in the fiducial profile and the 6.68 au gap needs d = 0.92 as opposed to d = 0.81.

With the inclusion of the second ring, the resulting model image matches the 2021 June image as well as the ASB (see Figures 8 and 9). We therefore demonstrate that two mutually inclined rings at different PAs well represent the observed disk images and represent the best explanation for the change in behavior. At earlier times, the PAs of the two rings were closer together and their resultant shadows overlapped, causing an apparent single shadow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.4. Observational Implications for Two Tilted Rings

We report new images of the TW Hya protoplanetary disk in visible total-intensity scattered light. We have shown that the disk's azimuthal brightness asymmetries at r > 40 au, previously interpreted as arising from a misaligned inner disk and casting a single shadow, now appears to be comprised of two shadows, implying that a single disk at r < 6 au is no longer a favored explanation as the originating shadowing structure. This is the first detection of a shadow splitting into more than one feature on a protoplanetary disk.

We instead consider two inclined rings at roughly the same orbital separation. This suggested explanation, when coupled with a series of gaps in dust density and moderate disk flaring, matches both the one-dimensional SED as well as the detailed features seen in scattered light. Nonetheless, the models presented here are not necessarily unique since it is beyond the scope of this paper to fully sample parameter space. It is true, however, that a detailed joint fit of the SED, resolved submillimeter continuum images, interferometry data, and scattered light in both total and polarized intensity represents a significant number of constraints on both the disk surface structure and the presence of interior tilted rings. The presence of shadows adds unique constraints to such modeling, hopefully in a way that provides a more unique model for TW Hya's disk.

The recent behavior of the shadow throws into question the shadow rotation period of 15.9 yr inferred by D17. This period was originally assumed to be a precession period of the inclined disk, which required a fairly massive companion to be consistent with observations. Additionally, the recent discovery of a similar shadow in submillimeter CO line emission (Teague et al. 2022) suggests that, instead of pure rigid rotation, the shadow(s) may oscillate in PA due to the mutual precession of the disks. Future observations will reveal whether the behavior of the shadow is truly periodic, and what that period relates to physically. If the shadow angular motions can be better disentangled, it could open up more powerful predictions for physical processes driving the shadow motions.

The presence of multiple tilted rings is highly suggestive of multiple planetary companions in the innermost regions of the TW Hya disk. While planets have been searched for in NIR data at around 6 au by VB17, no candidates were found down to a few tens of Earth masses in the existing disk gaps, assuming the level of disk extinction and planet evolutionary state is well known. Additional unsuccessful planet searches have been conducted with the Keck Vector Vortex coronagraph (Ruane et al. 2017) and in Hα (Cugno et al. 2019; Follette et al. 2022). Intriguing compact sources in the submillimeter have also been pointed to as evidence for planetary companions in the outer disk (e.g., Tsukagoshi et al. 2019). Additionally, the limited radial velocity data taken for TW Hya in the NIR is not sufficient to test the prediction of giant planet mass companions at distances of 5–10 au by itself, and visible radial velocities are plagued by the significant accretion and stellar activity of TW Hya (Huélamo et al. 2008). The Gaia mission's DR3 astrometric residual noise is not unusual compared to other stars of similar magnitude. Limits from a combination of Hipparcos and Gaia EDR3 proper motion anomaly measures show that the motion is consistent with no observable acceleration and place limits of 1.32, 1.29, 1.30, and 6.69 M_Jup for companions at orbits of 3, 5, 10, and 30 au, suggesting that if any planets are causing the structures that they are likely less than these masses (Brandt 2021; Kervella et al. 2022). Given the inferred orbital separation of the shadowing structures, the favorable face-on nature of the disk, and the sensitivity that the full Gaia astrometric observations will have, it is possible an astrometric signature of the perturbing companions could be detected in the future. The combination of a measured planet mass and resulting shadowing structure would mean that TW Hya will remain a key object for the study of planet formation and disk–planet interactions.

Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via 10.17909/mtx2-vz53. We thank the anonymous referee for their thoughtful comments. We additionally thank Emily Rickman for helpful discussions related to companion limits to TW Hya and Inbok Yeah for discussions related to MCFOST modeling. Support for program #16228 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This work is also 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 programs #11608, 13753, and 13665. The authors wish to thank the STScI Program Coordinator, S. Meyett, for her tireless work in getting the visits for Program 16228 scheduled successfully. R.N. acknowledges funding from UKRI/EPSRC through a Stephen Hawking Fellowship (EP/T017287/1). R.A. gratefully acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 681601), and from the Science & Technology Facilities Council (STFC) through Consolidated grant No. ST/W000857/1. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.