Chasing Shadows: Rotation of the Azimuthal Asymmetry in the TW Hya Disk^*

TW Hya is a nearby protoplanetary disk-hosting star that has been heavily studied over a broad range of wavelengths (e.g., Brickhouse et al. 2010; France et al. 2012; Kastner et al. 2014; Nomura et al. 2016). Its nearly pole-on disk (i < 7°; Qi et al. 2004) displays indications of active planet formation, despite its inferred age of 7–10 Myr (Weinberger et al. 2013; Ducourant et al. 2014; Herczeg & Hillenbrand 2014). Features in the disk that are indicative of such planet formation include an inner clearing in the disk out to ∼2.6 au for millimeter-sized dust (Andrews et al. 2016), a partially cleared gap at 89 au in scattered optical and near-IR (NIR) light (Krist et al. 2000; Weinberger et al. 2002; Debes et al. 2013), and a gap at 22–27 au observed in polarized scattered light (Akiyama et al. 2015; Rapson et al. 2015). A new interior gap has also recently been reported at 7 au by van Boekel et al. (2016).

We note that TW Hya's distance as determined by Gaia, 59.5 pc (Gaia Collaboration et al. 2016a), is roughly 10% larger than, but within 1σ of, that determined from Hipparcos data by van Leeuwen (2007) using the Hipparcos mission. The distance for TW Hya is still within 1σ of the Hipparcos value. We have thus adopted the Gaia value for the parallax, and all linear distances in au reflect this slight change compared to previous discussions in the literature.

We have reobserved the TW Hya disk using the newly commissioned BAR5 occulter on the Hubble Space Telescope/Space Telescope Imaging Spectrograph (HST/STIS) to achieve an inner working angle (IWA) of 02, which rivals that of ground-based polarized light observations. The images push to smaller angular distances from the host star than previous STIS images of the disk taken at the WEDGEA1.0 position (Roberge et al. 2005), enabling the exploration of the inner disk structure in scattered light. In addition to pushing inward in visible light, these new STIS images allow us to investigate whether there are any time-variable features in the disk's surface brightness distribution, as has been observed in other nearby disks (e.g., MP Mus; Schneider et al. 2014; Wolff et al. 2016).

In Section 2 we describe the observations of TW Hya with STIS, while in Section 3 we present an investigation of the inner surface brightness (SB) of the disk and report the apparent rotation of an azimuthal asymmetry present in the disk. In Section 4 we determine some possible physical mechanisms that could account for the variability we observe, and in Section 5 we discuss the implications of our observations.

Two new images of TW Hya were obtained with the STIS CCD on 2015 May 24 and 2016 March 24 (1024 × 1024 pixels, 50.7 mas pixel^–1), separated by roughly a year (GO #13753; PI: Debes). Table 1 lists all TW Hya observations used in this manuscript. In each new STIS epoch, two contiguous orbits targeting TW Hya were executed, immediately followed by one orbit targeting a point-spread function (PSF) reference, HD 88512, and finally an orbit of TW Hya. We used the same PSF star with earlier epochs of STIS data, to remove one source of uncertainty if temporal changes were observed. TW Hya was observed at three different sky orientations separated by ±15° at each epoch.

Table 1.  New and Archival HST Observations of TW Hya

MJD	UTC Date	Instrument	Filter	Ref. Star	F${}_{\star }$/F${}_{\mathrm{ref}}$	T${}_{\exp }$ (s)	Orientatione
51,041	1998 Aug 16	NICMOS	F160W	GL 879	0.031	1213	−72.7, −65.7
51,142	1998 Nov 25	NICMOS	F110W	${\tau }^{1}$ Eridani	0.099	1213	52.5, 59.5
51,851	2000 Nov 03	STIS	50CCD	HD 85512	0.055	2164	−139.9
51,905	2000 Dec 27	STIS	50CCD	HD 85512	0.057	2164	−106.9
53,037	2004 Feb 02	NICMOS	POL*La	BD+32°3739	3.5b	863	88.6, 118.5
53,108	2004 Apr 04	NICMOS	POL*L	BD+32°3739	3.4	863	182.7, 212.6
53,538	2005 Jun 17	NICMOS	...c	CD-43 2742	...d	3712	−125.4, −97.4
57,166	2015 May 24	STIS	50CCD	HD 85512	0.051	5661	49.0, 64.0, 79.0
57,471	2016 Mar 24	STIS	50CCD	HD 85512	0.055	6381	−31.9, −16.9, 3.8

Notes.

^aThe POL filters are POL0L, POL120L, and POL240L. ^bThe scaling factors for the POL filters were averaged over all spacecraft orientations and filters. ^cF171M, F180M, F222M. ^dSee Debes et al. (2013) for details of the scaling for these filters. ^eAll values in this column are in degrees.

Download table as:  ASCIITypeset image

Within each orbit, the star was placed behind the BAR5⁷ (width = 015) occulter. A three-point dither was executed with a width of ±0.25 pixels (±13 mas) in the direction perpendicular to the long axis of the occulter, and then the star was placed behind the WEDGEA1.0 (width = 05) occulting position. Dithering mitigates contrast loss due to target acquisition nonrepeatabilities. Exposures at BAR5 were short, to ensure that the star did not exceed the CCD full well. We utilized subarrays to minimize overheads due to readouts, with the BAR5 subarray consisting of a 1024 × 100 (519 × 5071) section of the detector centered on TW Hya and HD 85512. The WEDGEA1.0 subarrays consisted of 1024 × 427 (5193 × 2165) sections of the detector. In the first epoch, we obtained 54 exposures with integration times of 17 s each, and nine exposures with 527s integration times at WEDGEA1.0. Because the light in the PSF at the edge of the bar was smaller than originally estimated by roughly a factor of 5, we increased the BAR5 exposure times for the second epoch, such that there were now 18 exposures of 91 s each. For reference, the average observed count rate per pixel at the edge of BAR5 on the unfiltered STIS CCD with GAIN = 4 for TW Hya was 140 counts s⁻¹. The WEDGEA1.0 exposure times were unchanged for the second epoch. Total exposure times for each epoch are listed in Table 1.

TW Hya was also observed over two epochs separated by 2 months in 2000 with STIS/WEDGEA1.0, also using HD 85512 (GO #8624; PI: Weinberger, Roberge et al. 2005). The spacecraft orientation difference between the two 2000 epochs was 33° for the science target. In 2000, a background star near HD 85512 is present at ∼08. At one orientation it is near a diffraction spike, whereas at the second orientation it is close to the inner region of the disk. By 2015–2016, the background star had moved >4'' away, due to HD 85512's high proper motion. We isolated the star from HD 85512's PSF by subtracting one 2000 epoch from the other and creating an empirical PSF from the combined image. We then subtracted this template from the original images. We also masked other background stars that imprinted negative residuals onto our PSF-subtracted images.

We aligned HD 85512's position to that of TW Hya and subtracted off a scaled reference PSF. We centroided images by linearly fitting the vertical position of the diffraction spikes and solving for the intercept of the spikes to an accuracy of better than 0.1 pixels (Schneider et al. 2014). Further refinements to the measured offsets and scalings were obtained by iteratively subtracting the scaled PSF from the target and minimizing a ${\chi }^{2}$ metric for centroid shifts and scalings. For BAR5, we minimized subtraction residuals at the mask edge by selecting the best matches in dither position. Finally, we applied a custom coronagraphic flat field to all data, which was constructed from available flat-field images of the 50CORON aperture.⁸ This corrects flux right near the occulters at the few percent level.

We reverified the intensity scaling relations used between HD 85512 and TW Hya in Roberge et al. (2005). In the STIS passband, TW Hya has been variable at the level of 5%–10%. The scalings used in PSF subtraction are presented in Table 1. We confirmed that TW Hya is not variable by greater than 2% orbit to orbit via target versus target subtractions in our 2015 and 2016 epochs. We applied the stellar intensity scaling determined from the WEDGEA1.0 observations also to those at the BAR5 position where the disk fills the subarray. For the 2015 and 2016 epochs, we combined the WEDGEA1.0 and BAR5 PSF-subtracted images by averaging overlapping areas near the occulters out to 2''. Missing data at one position were filled with data from the other, if any existed.

STIS high-contrast images can be particularly sensitive to color mismatches, which tend to impact disk photometry with azimuthally symmetric zones of over- or undersubtraction of the target PSF interior to ∼1'' (Grady et al. 2003; Schneider et al. 2014). HD 85512's spectral energy distribution does not exactly match TW Hya across the STIS CCD bandpass, but the impact to PSF subtraction appears negligible at the level of <10%, given the good agreement in the slope of the surface brightness at these distances compared to those measured in the NIR with NICMOS, GPI, and HiCIAO (see Section 3).

We noted that there existed a small difference in fractional disk brightness between BAR5 and WEDGEA1.0. This is likely due to differing charge transfer efficiency at each occulter location from radiation damage on the STIS/CCD detector (Anderson & Bedin 2010). The exact impact of these effects is dependent on the clocking of the detector, the temperature of the detector, the amount of charge within each pixel, the distance each row of the CCD is from the amplifier, and the background present. Given the complexity of such an issue, we did not try to quantify these effects. Mitigating this problem by scaling up the WEDGEA1.0 position images by 5% allowed for a better empirical match between the BAR5 and WEDGEA1.0 images. The effect adds an overall 5% systematic uncertainty to the absolute SB of the disk at radii where BAR5 and WEDGEA1.0 overlap.

In addition to our STIS observations, we also used NICMOS images from programs 7226 (PI: Becklin), 10167 (Pi: Weinberger), and 9786 (PI: Hines). The last program in part consists of 2 $\mu {\rm{m}}$ coronagraphic imaging polarimetry of TW Hya from 2004. Details describing the NICMOS polarimetric observations of TW Hya, as well as an in-depth analysis of the Stokes images, will be presented elsewhere (C. A. Poteet et al. 2016, in preparation). In this study, however, we utilize the 2 $\mu {\rm{m}}$ polarimetric total intensity image (Stokes I) of TW Hya and briefly describe its construction.

Following the standard observing strategies for NICMOS coronagraphy, TW Hya was observed in 2004 February and April, and the unpolarized PSF reference star BD +32°3739 was observed between 2003 November and 2004 August with the NICMOS Camera 2 linear polarizers (POL0L, POL120L, and POL240L). The data were processed using high-level science products from the Legacy Archive PSF Library And Circumstellar Environments database (LAPLACE; Schneider et al. 2010), which are recalibrated using improved techniques in flat-field correction, dark subtraction, and bad pixel correction. During our processing procedure, the images were optimally registered using only the diffraction pattern produced by the HST secondary mirror support. In addition, the data were co-aligned in the science instrument aperture frame, background-subtracted, and further examined for the presence of bad pixels. Finally, the visit-specific co-aligned images were median-combined for each linear polarizer.

To remove the PSF component from the TW Hya polarimetric data, classical PSF subtraction was implemented using the unpolarized star BD +32°3739. After masking the central region surrounding the coronagraph, the scale factor between the TW Hya and BD +32°3739 median-combined images was determined for each linear polarizer by minimizing the residual signal along the full extent of the diffraction pattern. Adopting the polarimetry reduction algorithm of Hines et al. (2000), co-oriented Stokes I, Q, and U images were produced from the visit-specific PSF-subtracted images and subsequently smoothed using a 3-pixel FWHM Gaussian. The resulting images for all visits were then median-combined to construct the final Stokes images.

Figure 1 shows logarithmically scaled images of the year 2000, 2004 (NICMOS POL*), 2015, and 2016 epochs. The last two epochs have full angular coverage down to ∼05 (29 au) from the star, and partial angular coverage down to ∼02 (11.8 au). For the STIS images we assumed a photometric conversion of 8.43e–7 Jy count⁻¹ s pixel⁻¹ at an effective wavelength of 7204 Å as given by a Pysynphot calculation (STScI Development Team 2013) of STIS's sensitivity for an M0 spectral type. We choose an M0 spectral type for TW Hya in this wavelength range because this most closely represents the combination of the stellar photosphere and its veiling due to accretion (Debes et al. 2013). In 2000 the photometric conversion is 1.6% lower, due to higher sensitivity of the instrument. We present images in Figure 2 where the brightness within each pixel is scaled by the square of the stellocentric distance to highlight features within the disk. The inner part of the disk appears dim in these images because the disk SB profile is shallow from the inner working angle of the images out to 07 (38 au). Between 07 and 1'' (42–59 au), a ring is present, with clear azimuthal asymmetries. A second radial gap is clearly seen at 15 (89 au), with a second fainter ring at 23 (137 au).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We measured azimuthally averaged SB profiles of TW Hya's disk from 2000, 2015, and 2016. The uncertainty in our PSF scaling translates into at most a 5% systematic uncertainty in disk SB per pixel. Any residual systematic uncertainties, which dominate over statistical uncertainties, were accounted for by measuring the standard deviation in our photometric apertures. We compare our SB profiles to those seen in polarized NIR light and to the ALMA dust continuum maps at high resolution. In addition to our STIS observations, we also use NICMOS images from programs 7226 (PI: Becklin) and 10167 (PI: Weinberger).

3.1. Radial Surface Brightness of TW Hya's Disk

The left panel of Figure 3 shows the 360° azimuthally averaged SB profiles for the 2000 epoch compared to a combination of 2015/2016 taken with STIS. We recover the broad behavior of the disk as seen in previous work (Akiyama et al. 2015; Rapson et al. 2015; Debes et al. 2016). In particular, we confirm the presence of a depletion of SB from our inner working angle of 10 au to 38 au. Following Akiyama et al. (2015), we report the STIS SB power-law fits to the region between 02 and 04 and the region between 04 and 08. We find slopes of −1.2 ± 0.2 and −0.4 ± 0.1, compared to −1.7 and −0.4 for the recent GPI images of the disk (Rapson et al. 2015) and −1.4 and −0.3 for the recent HiCIAO image of the disk (Akiyama et al. 2015). The right panel of Figure 3 shows the excellent match between the total intensity light in the visible and the polarized intensity light at J and K.

Zoom In Zoom Out Reset image size

We averaged the images from the 2015 and 2016 epochs to determine the extent of the disk (left panel of Figure 3). The resulting azimuthally averaged profile shows disk flux at 10σ out to 76, or 452 au, which is larger than previously reported (Roberge et al. 2005). Beyond this point, the disk is not detected significantly, due to the systematic uncertainty in the background level of the images, which is denoted by the dashed lines in the left panel of Figure 3. Using our above flux conversions, the SB of the disk edge is 0.2 μJy arcsec⁻² at 76. Comparing this limiting SB to the star, the per-pixel contrast of the disk relative to the flux density of the star (0.18 Jy; Debes et al. 2013) at this angular distance is 3 × 10⁻⁹.

Over 15 yr, the SB profile of the disk may have varied, especially if there were any systematic changes in the illumination due to self-shadowing from the inner part of the disk or from the central star itself. We investigated whether there were any changes between the 2000 and 2015/2016 epochs (left panel of Figure 3). The disk SB does not appear to vary by more than 5%–10% between epochs from the inner working angle to the outer edge of the disk.

In addition to variability, we looked for radially coincident features in the TW Hya disk scattered-light emission and the ALMA brightness temperature (right panel of Figure 3; Rapson et al. 2015; Andrews et al. 2016). The breaks seen in SB are coincident with the 27 au gap reported in the ALMA data, which is just outside the submillimeter transition between the optically thick inner disk and the optically thin outer disk (Andrews et al. 2016; Nomura et al. 2016). The gap might be caused by a planet, especially since there is a lack of millimeter-sized dust observed at 22 au (Tsukagoshi et al. 2016). Alternatively, freezing out and grain growth can cause the apparent gap. As the opacity to stellar irradiation decreases, this causes shadowing at the disk surface and cooling of the disk locally (i.e., Andrews et al. 2016; Debes et al. 2016). Since stellar irradiation is the primary heat source in the disk at this distance, the reduction in temperature can propagate down to the midplane of the disk, especially if the disk is optically thin (Jang-Condell & Turner 2012). Since the submillimeter continuum observations are sensitive to both disk temperature and density, such a temperature reduction would lower the local surface brightness. The gap is likely to be narrower and deeper in the case of a density decrement as opposed to a pure temperature decrease, but fitting an exact model is beyond the scope of this paper. Apart from the region near the 24 au gap, there are no clear correlations between the STIS SB and the submillimeter emission.

3.2. Azimuthal Asymmetry Changes

The TW Hya disk has a sinusoidal azimuthal asymmetry at certain radii (Roberge et al. 2005), with a peak position angle that does not match the position angle of the semimajor disk axis of 155° (Rosenfeld et al. 2012; Andrews et al. 2016). The asymmetry was further observed at larger radii in the same STIS data set and also detected in an earlier NICMOS F110W image (Debes et al. 2013). The origin of the SB asymmetry is debated and has been attributed to inclination effects for flared disks along the disk minor axis (Debes et al. 2013), the forward scattering properties of the disk (Roberge et al. 2005), or the presence of a warp in the disk that modified the azimuthal surface brightness profile (Roberge et al. 2005; Rosenfeld et al. 2012).

If any inner disk structure blocks the light of the star, this can cause shadowing of the outer disk that modifies the azimuthal SB profile of a disk. An inner warp may also be a signature of a planet, since flaring and temperature variations in the disk are expected at the edges of planet gaps (Jang-Condell & Turner 2012). In this section, we demonstrate that the position angle of the SB asymmetry peak in the TW Hya disk is not stationary, but moves with time at disk radii >50 au, placing constraints on the origin of the azimuthal scattered-light asymmetry. Further, we show that interior to 50 au, there is a change in the nature of the SB azimuthal asymmetry.

To measure the properties of the asymmetries, we divided each azimuthal SB profile for a given radius by the mean SB at that radius to allow for unbiased radial averaging. We chose radii that matched seven regions of the disk centered at 046 (27 au), 066 (39 au), 089 (53 au), 114 (68 au), 147 (88 au), 183 (109 au), and 236 (141 au), with widths of 02 (12 au), 02 (12 au), 025 (15 au), 025 (15 au), 041 (24 au), 030 (18 au), and 056 (33 au) respectively. These were chosen to roughly correspond to regions within and near the two gaps seen in scattered light at 22 and 88 au. Figure 4 shows a comparison between our 2015 and 2016 data at the above distances, where azimuthal asymmetries are clearly detected in the data at all radii. We first tested each radius to see whether it had a significant departure from a constant profile by calculating its ${\chi }_{\nu }^{2}$ value. If the ${\chi }_{\nu }^{2}$ value of the azimuthal profile had a less than 10⁻³ probability of being consistent with a flat profile, we fit a cosine function of the form

$$\begin{eqnarray}&&F(\theta )=A\cos [2\pi (\theta -{\theta }_{{\rm{o}}})]+B,\end{eqnarray} \tag{ 1 }$$

where θ is the position angle as measured east through north on the sky, A is the amplitude of the asymmetry, ${\theta }_{{\rm{o}}}$ is the position angle of the asymmetry peak, and B is the constant offset to the curve.

Zoom In Zoom Out Reset image size

In all cases the resulting fit had an improved ${\chi }_{\nu }^{2}$ value that was close to 1. Table 2 shows our best-fit cosine parameters for each epoch and radius measured.

Table 2.  Summary of Best-fit Asymmetry Parameters

MJD	Filter	R (au)	Amplitude	PA (deg)	Offset	${\chi }_{\nu ,\mathrm{fit}}^{2}$	ν
51,142	F110W	38	0.35 ± 0.12	192 ± 4	1.15 ± 0.08	0.62	5
51,142	F110W	51	0.24 ± 0.03	206 ± 4	0.99 ± 0.02	1.87	9
51,142	F110W	67	0.25 ± 0.02	207 ± 4	1.00 ± 0.02	1.33	10
51,142	F110W	87	0.18 ± 0.02	206 ± 5	1.01 ± 0.01	3.02	13
51,142	F110W	112	0.20 ± 0.02	207 ± 3	1.02 ± 0.01	5.73	15
51,142	F110W	141	0.21 ± 0.01	211 ± 3	1.03 ± 0.01	13.02	15
51,878	STIS	39	0.20 ± 0.02	243 ± 3	1.02 ± 0.01	3.04	7
51,878	STIS	53	0.20 ± 0.02	256 ± 3	0.97 ± 0.01	0.77	7
51,878	STIS	68	0.21 ± 0.01	271 ± 3	0.98 ± 0.01	1.91	10
51,878	STIS	88	0.28 ± 0.01	257 ± 2	0.97 ± 0.01	2.35	15
51,878	STIS	109	0.21 ± 0.02	253 ± 3	1.02 ± 0.01	0.48	15
51,878	STIS	141	0.07 ± 0.01	223 ± 7	1.00 ± 0.01	1.84	15
53,072	POL*	38	0.14 ± 0.02	260 ± 7	1.00 ± 0.01	1.55	15
53,072	POL*	51	0.20 ± 0.02	279 ± 5	1.03 ± 0.01	0.98	15
53,072	POL*	67	0.34 ± 0.03	284 ± 3	1.01 ± 0.02	1.95	15
53,072	POL*	87	0.21 ± 0.03	296 ± 6	1.01 ± 0.02	0.85	15
53,072	POL*	112	0.13 ± 0.02	314 ± 11	1.02 ± 0.02	1.78	15
53,072	POL*	141	0.14 ± 0.04	16 ± 19	0.99 ± 0.02	2.00	15
53,538	Various	51	0.09 ± 0.02	353 ± 12	1.01 ± 0.01	2.72	15
57,166	STIS	27	0.25 ± 0.03	228 ± 7	1.01 ± 0.02	0.71	14
57,166	STIS	39	0.19 ± 0.01	232 ± 2	1.00 ± 0.01	0.89	15
57,166	STIS	53	0.13 ± 0.01	215 ± 3	0.99 ± 0.01	1.82	15
57,166	STIS	68	0.13 ± 0.01	219 ± 4	0.99 ± 0.01	1.53	15
57,166	STIS	88	0.21 ± 0.01	229 ± 3	0.98 ± 0.01	1.33	15
57,166	STIS	109	0.18 ± 0.01	229 ± 2	1.01 ± 0.01	0.78	15
57,166	STIS	141	0.10 ± 0.01	222 ± 4	1.01 ± 0.01	0.65	15
57,471	STIS	27	0.18 ± 0.03	220 ± 8	1.02 ± 0.02	0.47	14
57,471	STIS	39	0.20 ± 0.01	233 ± 2	0.99 ± 0.01	0.52	15
57,471	STIS	53	0.14 ± 0.01	234 ± 3	0.99 ± 0.01	1.07	15
57,471	STIS	68	0.10 ± 0.01	229 ± 4	1.00 ± 0.01	2.21	15
57,471	STIS	88	0.18 ± 0.01	243 ± 2	1.02 ± 0.01	1.05	15
57,471	STIS	109	0.17 ± 0.01	245 ± 1	1.01 ± 0.01	2.43	15
57,471	STIS	141	0.12 ± 0.01	242 ± 3	1.01 ± 0.01	0.76	15

Download table as:  ASCIITypeset image

Figure 5 shows that between 2015 and 2016, the position angle of the azimuthal brightness peak shifted by a significant amount, on average 17° ± 4°, beyond 50 au. Interior to 50 au, there is no evidence of a change in the position angle of the asymmetry peak. Since the two epochs are roughly separated by 10 months, this translates to an instantaneous angular velocity of 20° ± 5° yr⁻¹.

Zoom In Zoom Out Reset image size

To further determine whether the azimuthal asymmetry was variable, we constructed similar SB asymmetry profiles from the 2000 STIS image, the 1998 F110W/F160W NICMOS images, the 2004 POL*L NICMOS total intensity (Stokes I) image, and the 2005 F171M/F180M/F222M NICMOS images. We analyzed the NICMOS images in a similar fashion to the STIS data.

For each filter we determined whether an asymmetry was significantly detected by calculating ${\chi }_{\nu }^{2}$ values as described for the two STIS epochs in 2015 and 2016. With the exception of the F160W image, we found significant asymmetries for at least one radius at every filter wavelength, implying that the asymmetry is present from the optical through the NIR. The F160W images show no appreciable asymmetry, and we conducted tests that show that the systematic noise within the image prevents the secure detection of any similar asymmetry seen at other wavelengths.

We also investigated whether a systematic PSF residual feature could masquerade as a moving asymmetry by measuring the position angle asymmetry in single spacecraft orientations. For all of the STIS data, the asymmetry is fixed in sky orientation at a given epoch, i.e., the orientation of the asymmetry is not correlated with changing spacecraft orientation. The two F110W spacecraft orientations were only separated by 7°, but the average of the position angles is consistent with the asymmetry being fixed on the sky. For the NICMOS medium band and POL data, the dependence (or lack thereof) on spacecraft roll is less well determined, due to lower signal-to-noise ratio. If there were systematic residual features in the medium band and POL data, they would have to reside at roughly 100° in the medium-band data and 140° in the POL data to mimic the behavior of the asymmetry's motion.

Figure 6 shows the azimuthal profiles for five of the epochs with significantly detected asymmetries at 88 and 109 au. The position angle of the peak in the azimuthal surface brightness profile clearly moves with time. Sampling of the asymmetry is nonexistent between 2005 and 2015, but the motion is consistent with constant rotation of the peak, as would be expected for circular orbital motion of a feature. In Figure 7 we plot the position angle of the asymmetry peak as a function of time at 53, 68, 88, 109, and 141 au. If we take an average of the position angles over radius and assume that the scatter in the data represents the true uncertainty in the position angle measurement, a linear least-squares fit to the data results in a rate of 227 ± 04 yr⁻¹ in a counterclockwise direction, consistent with periods of 15.9 ± 0.3 yr. The measured velocities greatly exceed the local Keplerian velocity at the radii sampled (53–141 au; expected angular velocities of ∼08−02 yr⁻¹), but would be consistent for an object orbiting at 5.6 au.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The recently detected asymmetries in SPHERE data look qualitatively different from our total intensity images (van Boekel et al. 2016). It is difficult to perform a direct comparison with these results since polarized intensity profiles may have different geometrical effects than those of total intensity images. For example, the 2015 data from SPHERE were taken within a few months of our STIS data. While the total intensity asymmetry does not appear to significantly change amplitude or shape with wavelength according to our data, the SPHERE data show a series of smaller peaks at roughly 130° and 220°, with a decrement of flux at 300°.

The angular velocity of the asymmetry is larger than the expected Keplerian velocity by at least an order of magnitude and shows a constant velocity across a large swath of radii between 50 and 141 au, yet no motion interior to 50 au. A precessing warped inner disk that is at a higher inclination than the outer disk and is shadowing parts of the outer disk could explain the observed behavior in TW Hya. A warped inner disk at ∼5 au was invoked to explain anomalous CO line emission in the disk, as well as the azimuthal asymmetry seen in the disk SB (Rosenfeld et al. 2012). The recent detection of a gap at ∼7 au in scattered light is close to this location (van Boekel et al. 2016). There has also been a tentative detection of a point source at  6 au with nonredundant aperture masking of the star (Willson et al. 2016). The period of the asymmetry motion, 15.9 yr, is comparable to the expected period from a circular orbit around TW Hya at 5.6 au.

We consider whether an inclined inner disk that shadows the outer part of the disk while precessing owing to a planetary companion is a plausible explanation for what we observe. The inclination of the inner disk may not be significantly different from that of the outer disk, contrary to what is seen in HD 142725 (Marino et al. 2015), which only shows sharp shadows along the line of nodes for the inner, highly inclined disk. Moderately inclined inner disks that were sufficiently flared or puffy could result in an asymmetry that roughly covers half of the outer disk, corresponding to where one lobe of the inner disk shadows the flared surface of the outer disk. This is what we see for TW Hya.

The ALMA 870 μm continuum emission image of TW Hya (Andrews et al. 2016) provides some possible constraints on interpreting the location of the shadowing structure. Since we need to invoke an inclined inner disk, we note that the most likely location is coincident with the unresolved emission detected around TW Hya interior to 1 au. This inner disk has also been inferred to exist in interferometric data sets in the NIR (Akeson et al. 2011). There are no obvious warps or departures from azimuthal symmetry in the submillimeter beyond 1 au, and the interior of the system is preferable for shadowing a large range of grazing angles in the outer disk.

To fully determine the plausibility of such an inclined disk and its impacts on the outer disk, detailed hydrodynamical modeling and radiative transfer would need to be done to self-consistently include the features observed in thermal emission and scattered light. While such modeling is beyond the scope of this paper, we can assume that the outer disk is well approximated from the best-fit model of TW Hya presented in Debes et al. (2013), and that the inner inclined disk has an opening angle that is constrained by the extent of the shadow asymmetry. The inner disk could be torqued by a companion on an inclined orbit, and it will precess owing to the interaction between the companion's orbit, the spin of the star, and the disk (e.g., Lai 2014). Based on the best-fit model of TW Hya's disk in Debes et al. (2013), the height H of the $\tau =2/3$ surface above the midplane for the disk in the STIS bandpass is 12 au at a radius R = 53 au and H = 34 au at R = 141 au, the rough extent to which we recover the shadow asymmetry in all three STIS epochs. We determine the range of inclinations relative to the disk midplane that cause the shadow by calculating $\delta ={\sin }^{-1}(H/R)$. If we assume that the inner disk midplane is in the middle of the angles probed, the inferred inclined disk must span from ∼13° to ∼14°. We note that this is compared to the plane of the outer disk–in reality the outer disk is inclined from our line of sight by 7° (Qi et al. 2004).

The lack of rotation in the asymmetry interior to 53 au could be due to several factors: (1) the opening angle of the inner disk is larger and the difference in inner versus outer disk inclinations is smaller, such that the inner disk casts shadows over all azimuthal angles; (2) the outer disk is fully illuminated by the star, and thus there exists a puffed inner rim at ∼2 au, shadowing the disk from 2 to 53 au; or (3) the inner edges of one of the other gaps (13 or 22 au) are shadowing the region between 27 and 53 au. In any case, the asymmetry present at these radii has a peak at a position angle of 232°, close to the expected minor axis of the disk. Given this alignment, the origin of the (stationary) asymmetry interior to  50 au presumably lies in a combination of inclination effects and the proximity to the gap at  25 au discussed previously (Akiyama et al. 2015; Rapson et al. 2015; Debes et al. 2016).

The period of rotation of the asymmetry, if caused by precession of an inclined inner disk due to a planetary companion, is dictated by the mass of the companion, the extent of the perturbed disk, and the orbit of the companion. Given the preceding model for the rotating outer disk shadow pattern that invokes an inclined inner disk, we can place tentative constraints on the mass of the putative perturbing body. The precession period is given by (Lai 2014)

$$\begin{eqnarray}&&P=\displaystyle \frac{2.04\times {10}^{5}}{{\mu }_{{\rm{c}}}\cos {i}_{{\rm{c}}}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{{r}_{\mathrm{disk}}}{1\mathrm{au}}\right)}^{-\tfrac{3}{2}}{\left(\displaystyle \frac{{a}_{{\rm{c}}}}{300\mathrm{au}}\right)}^{3}\,\mathrm{yr},\end{eqnarray} \tag{ 2 }$$

where ${\mu }_{{\rm{c}}}$ is the mass ratio of the perturber to the central star, i_c is the inclination of the binary relative to the inner disk, ${M}_{\star }$ is the mass of the central star (0.69 ${M}_{\odot }$), r_disk is the extent of the inner disk, and ${a}_{{\rm{c}}}$ is the semimajor axis of the perturber. Solving for the perturber mass with i_c = 14°, P = 15.9 yr, ${r}_{\mathrm{disk}}=1$ au, and a_c = 1.1 au results in a mass of 1.2 M${}_{\mathrm{Jup}}$. If there were a companion at a_c = 7 au, with ${r}_{\mathrm{disk}}=6$ au, the mass of the companion would need to be 19 M${}_{\mathrm{Jup}}$ to create a 15.9 yr precession period in the inner disk. The high mass needed at this distance suggests that the location of the inclined disk is closer to 1 au than at this location if disk precession is the cause of the rotating shadow.

The radial velocity semiamplitude and period for such a planetary companion at 1.1 au would be ≈10 m s⁻¹ and ≈1.4 yr, respectively, which is consistent with the nondetection of any radial velocity trend in TW Hya at 6 m s⁻¹ precision over a 6-day period (Huélamo et al. 2008; Figueira et al. 2010). The proposed structure of the disk here would also need to be reconciled with existing interferometric data in the NIR, in the mid-IR, and at centimeter wavelengths, as well as with the reported gap in scattered light reported by van Boekel et al. (2016).

The sparse temporal coverage of the extant HST scattered-light imaging of the TW Hya disk leaves room for other potential behavior for the asymmetry that may mimic rotation at constant angular velocity. It is possible that the inferred shadow pattern rotation is due either to a much shorter period phenomenon or to some stochastic variation. However, we have inspected the individual STIS images obtained in 2000 and the individual STIS images of 2015 and 2016, and we find no significant variations on such (month to year) timescales. If we assume that we would have been able to detect azimuthal brightness pattern motion in these images at the level of 5°–10° per epoch, then the lack of detectable changes within individual epochs appears to rule out the subset of stochastic or very rapid variations that might occur on hourly to monthly timescales.

We can also rule out the possibility that the observed linear trend in position angle with time is the result of observations at a completely random set of position angles. We constructed a Monte Carlo model of random position angles as a function of time at the six epochs in which we observed TW Hya. We performed 10⁶ trial slope fits by taking randomly selected position angles as a function of the epochs we sampled. We determined what fraction of these trials could successfully be fit with a line with the criteria that a X² probability of success was better than 5%. From these 10⁶ trials, we find a 4 × 10⁻⁵ probability that the slope is a random occurrence.

While it is possible that the asymmetry is only due to a shadow, inclination effects can also modify the inferred motion of a shadow asymmetry. As mentioned in Section 3.2, the inclination of a flared disk can result in a weak azimuthal asymmetry that can become more pronounced as the inclination departs from face-on or in the presence of forward-scattering dust grains. The combination of the shadow and the inclination asymmetries for small inclinations can be approximated analytically:

$$\begin{eqnarray}&&{{SB}}_{\mathrm{disk}}\approx {{SB}}_{{\rm{o}}}[1+{A}_{1}\sin \theta +{A}_{2}\sin (\theta +\alpha )],\end{eqnarray} \tag{ 3 }$$

where A₁ is the amplitude of the inclination asymmetry, A₂ is the shadow asymmetry amplitude, θ is the sky position angle, and α is the position angle as a function of time for the shadow. This simplifies to

$$\begin{eqnarray}&&{{SB}}_{\mathrm{disk}}\approx {{SB}}_{{\rm{o}}}[1+C\sin (\theta +\phi )],\end{eqnarray} \tag{ 4 }$$

where the amplitude of the asymmetry is

$$\begin{eqnarray}&&C=\sqrt{{A}_{1}^{2}+{A}_{2}^{2}+2{A}_{1}{A}_{2}\cos \alpha }\end{eqnarray} \tag{ 5 }$$

and the observed position angle of the asymmetry is given by

$$\begin{eqnarray}&&\phi ={\tan }^{-1}\displaystyle \frac{{A}_{2}\sin \alpha }{{A}_{1}+{A}_{2}\cos \alpha }.\end{eqnarray} \tag{ 6 }$$

If ${A}_{1}\gt {A}_{2}$, then one should expect the observed asymmetry to change in amplitude and even disappear at certain epochs, while librating about the observed position angle of the minor axis of the disk. Conversely, if ${A}_{1}\lt {A}_{2}$, the asymmetry will be constant with time and circulate at the angular velocity of the shadow. In the case of TW Hya, given the low expected inclination of the outer disk, A₁ is a factor of ∼3 smaller than the amplitudes actually measured. No obvious change in amplitude is observed epoch to epoch between the STIS images. There are also no significant deviations from a purely linear trend in the measured position angle of the asymmetry, as would be expected for a shadow amplitude that is slightly larger than the inclination asymmetry.

With STIS+NICMOS, we can probe the changing illumination of the protoplanetary disk of TW Hya with unprecedented detail on timescales of 17 yr. The new BAR5 occulter on STIS probes stellocentric distances comparable to the orbital radius of Saturn. While the SB profile of the disk has remained constant over 15 yr, the azimuthal asymmetry previously seen in the disk is variable in position angle and best described by the orbital motion of a shadowing structure in the inner parts of the disk.

These observations are part of a growing list of scattered-light variability observed in protoplanetary transition disks with rings and gaps. For example, PDS 66 showed variability in its SB in the optical (Schneider et al. 2014) on monthly timescales and hints at azimuthal movement of features with GPI polarized intensity data (Wolff et al. 2016). The transition disk J160421.7-213028, observed with SPHERE/ZIMPOL in the visible and with HiCIAO in polarized H-band intensity, also showed a localized dip in the SB of its disk that appeared to move by ∼12° yr⁻¹; however, the angular velocity of the shadow pattern is indeterminate, given that data were obtained at only two epochs (Mayama et al. 2012; Pinilla et al. 2015). Additional claims of variable shadowing have been put forth for HD 163296 (Wisniewski et al. 2008), but again the epoch-to-epoch variation was hard to characterize owing to a dearth of coverage.

The combination of the broad azimuthal extent of the asymmetry, its rotational motion, and the range of radii over which it appears is suggestive of the presence of an anomalously inclined structure in the disk interior to 5 au that shadows the outer regions of the disk. Given the inferred rotational timescale of the azimuthal asymmetry, we propose that it is the result of shadowing by an inclined disk of radius <1 au that is precessing with an orbital period of 15.9 yr. Adopting reasonable assumptions concerning the inclination of the disk and its precession due to an external perturber at 1.1 au, we estimate the mass of a putative perturbing body to be  ∼${M}_{\mathrm{Jup}}$. Such a mass is small enough to have thus far escaped detection by previous campaigns to measure radial velocity variations in TW Hya.

Our results for TW Hya demonstrate that the detailed structure of protoplanetary disks within 20 au impacts the scattered-light emission observed in the outer disk. Time-domain high-contrast scattered-light imaging with high signal-to-noise ratio and small inner working angles, such as the observations of TW Hya with STIS, GPI, and SPHERE, provides an approach complementary to that of future observations of disks that will be obtained with ALMA and the James Webb Space Telescope.

Support for this work was provided by NASA through grants HST-GO-13753 and by the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. The HST data used in this article are available at the MAST archive (at 10.17909/T9XG63). J.H.K.'s research on protoplanetary disks orbiting nearby young stars is supported by NASA Exoplanets Program grant NNX16AB43G to Rochester Institute of Technology. This research has made use of the VizieR catalog access tool and the SIMBAD database, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS, 143, 23. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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. This manuscript benefited from discussions with Steven Lubow.