Multi-epoch Direct Imaging and Time-variable Scattered Light Morphology of the HD 163296 Protoplanetary Disk

We obtained high-contrast H-band imaging of HD 163296 using the High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) instrument (Hodapp et al. 2008) along with the AO-188 system (Hayano et al. 2008, 2010) at the Subaru Observatory on 2011 August 3. We used a circular occulting mask having a diameter of 03, and observed the system in standard Polarized Differential Imaging mode at four wave-plate positions (0°, 225, 45°, 675). We obtained 72 frames using 30 s exposures, yielding a total of 18 complete wave-plate sets. We determined that 8 wave-plate sets had lower AO performance, and discarded them during the reduction of the data. We also obtained a short, direct H-band photometric observation of HD 163296, and determined that the source's brightness at this epoch was 5.62 ± 0.05 mag.

We reduced our observations using standard double differencing techniques, as described in Hashimoto et al. (2011). To briefly summarize, the two subimages of each frame contain an ordinary and an extra-ordinary image, which can be summed and subtracted from their 90° counterparts to create stokes parameter −Q, +Q, −U, and +U images. The Q and U frames were then rotated into a common orientation, corrected for instrumental polarization, and summed to create final Q and U images. We corrected these data for the presence of a residual polarized halo having the properties of p = 1.00 ± 0.05% and θ = 425 ± 15. Final polarized intensity (PI) imagery was created from the total Q and U data, using $\mathrm{PI}=\sqrt{{Q}^{2}+{U}^{2}}$, as shown in Figure 1.

Zoom In Zoom Out Reset image size

To further simplify the analysis of our imagery, we adopt the now common practice of assuming single scattering, and rotated all of the light that is polarized perpendicular to the star by the angle ϕ into a Q_ϕ image and all of the light that is polarized parallel to the star into a U_ϕ image as defined below (Schmid et al. 2006).

$$\begin{eqnarray}&&{Q}_{\phi }=Q\times \cos 2\phi +U\times \sin 2\phi \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{U}_{\phi }=Q\times \sin 2\phi +U\times \cos 2\phi .\end{eqnarray} \tag{ 2 }$$

The final Q_ϕ and U_ϕ imagery for HD 163296 are shown in Figure 1. Little coherent signal appears present in the U_ϕ image, which helps confirm that little residual instrumental contaminants remain in these data. Next, we computed a signal-to-noise (SN) image, following the procedure outlined by Ohta et al. (2016). In summary, we computed the noise by measuring the standard deviation of every pixel in each of the Q and U frames used to construct the final imagery, then divided by the square root of the number of frames. The resultant SN image is shown in Figure 1.

We also observed HD 163296 multiple times with several near-IR instruments on NASA's Infrared Telescope Facility (IRTF) and at Apache Point Observatory (APO). We observed HD 163296 using the SpeX spectrograph (Rayner et al. 2003) at IRTF in its short-wavelength mode (0.8–2.4 μm) and long-wavelength mode (2.3–5.5 μm) on 2011 July 31, 2016 May 4, and 2018 June 24. These observations are contemporaneous with the HiCIAO 2011 observation (Section 2.1), the Gemini/GPI observation (Section 5.1), and the second SCExAO/CHARIS observations (Section 2.3) respectively. We observed HD 163296 using the TripleSpec spectrograph (Wilson et al. 2004) at the APO 3.5 m telescope, covering a spectral range of 0.95–2.46 μm, on 2018 May 16. This observation is contemporaneous with the first SCExAO/CHARIS observation (Section 2.3). We observed the nearby A0V star HD 163336 to perform telluric corrections for both the SpeX and TripleSpec observations. These data were reduced and calibrated using the standard reduction packages Spextool and Triplespectool (Vacca et al. 2003; Cushing et al. 2004). We also observed HD 163296 with The Aerospace Corporation's BASS, which covers two wavelength bands from 2.9–6 μm and 6–13.5 μm respectively, on 2011 August 1. HD 163336 was observed with BASS to flux calibrate these data. The instrument and data reduction method are fully described in Wagner et al. (2015). These SpeX, TripleSpec, and BASS spectra are plotted in Figure 2.

Zoom In Zoom Out Reset image size

We observed HD 163296 on 2018 May 22 and July 1 at the Subaru Observatory with SCExAO coupled with the CHARIS integral field spectrograph operating in low-resolution (R ∼ 20) and broadband (1.13–2.39 μm) mode, covering the JHK filters simultaneously (Groff et al. 2015). For the May observations, the conditions were stable with 04 seeing and 6–7 m s⁻¹ winds. Our observations consisted of co-added 60.4 s frames totaling ∼30 minutes of integration time and covering a modest parallactic angle rotation (ΔPA = 14.8^o). Due to highly variable conditions for the July observations, we obtained shorter exposures (30.9 s) and removed roughly 50% of the frames with poor AO correction, yielding ∼40 minutes of data covering 309 of parallactic angle motion.⁴¹

We followed the standard setup used for SCExAO/CHARIS broadband observations (Currie et al. 2018a; Goebel et al. 2018a), using the Lyot coronagraph with the 217 mas occulting spot and bracketing our coronagraphic sequence with blank sky frames to remove sky emission and instrumental artifacts. We used satellite spots produced from a 25 nm modulation on SCExAO's deformable mirror for spectrophotometric calibration and image registration (Jovanovic et al. 2015). For data cube extraction, we utilized the least-squares algorithm from the CHARIS Data Reduction Pipeline (Brandt et al. 2017). Basic data processing, including sky subtraction, image registration, etc., follows methods used for recent SCExAO/CHARIS broadband studies (Currie et al. 2018a, 2018b; Goebel et al. 2018b).

Spectrophotometrically calibrating CHARIS data for pre-transitional disk sources like HD 163296 requires either observations of a separate spectral standard or contemporaneous near-IR spectra. We opt for the latter, using the IRTF/SpeX and APO/Triplespec data previously discussed in Section 2.2. The spectra show only minor differences between epochs.

We explored a range of point-spread function (PSF) subtraction approaches leveraging on angular differential imaging (ADI; Marois et al. 2006), spectral differential imaging (SDI; Sparks & Ford 2002), and combinations of the two (ASDI, e.g., Marois et al. 2014). We further considered a variety of PSF subtraction algorithms, including A-LOCI (Currie et al. 2012, 2018a), Karhunen–Loéve Image Projection (KLIP; Soummer et al. 2012), and classical PSF subtraction (Marois et al. 2006). The approach implemented for κ And in Currie et al. (2018a), using A-LOCI to subtract the PSF in ADI and then again to remove residuals in SDI mode, yielded the best speckle suppression while preserving the signal from the disk. Due to the limited parallactic angle motion of both data sets (especially in May) and the presence of the disk, we utilized large optimization zones for the ADI step, employed local masking in the SDI step, and imposed a rotation/magnification criterion of δ = 0.5–1.0 PSF footprints in both steps to construct a reference PSF (see Lafreniére et al. 2007). For both steps, we used a singular value decomposition (SVD) cutoff of 10⁻⁶ to solve the set of linear equations that result in the weighted reference PSF for each region of each data cube slice (see Currie et al. 2015).

Figure 3 shows broadband (wavelength-collapsed) CHARIS images of HD 163296 from SCExAO/CHARIS for the May (left) and July (right) epochs after removing the stellar PSF through both ADI and SDI. Despite poor field rotation (May data) or variable conditions (July data), we clearly detect the outer ring of emission seen in polarimetry, which appears as a sharply defined crescent defining the forward-scattering edge of the structure. Self-subtraction footprints due to both ADI and SDI flank the ring. In individual passbands, the disk is just marginally visible in J-band but is well separated from residual speckle noise in H and K.

Zoom In Zoom Out Reset image size

We defined a conservative lower limit to the signal-to-noise ratio (S/N) of the trace of the disk in broad band, adopting the standard practice of replacing each pixel by the sum within its aperture, defining a radial-dependent noise profile, and applying a finite-element correction for the noise (Currie et al. 2011; Mawet et al. 2014). To be conservative, we include signal from the disk in our estimate of the noise profile. Except at the semiminor axis, where the disk signal is attenuated by self-subtraction, the disk trace is decisively detected, with a S/N per resolution element ranging from 3 to 8.5.

Our data do not reveal the candidate protoplanet identified in Keck/NIRC2 L_p data from Guidi et al. (2018) nor the companions predicted from ALMA data (Teague et al. 2018). The inner disk seen by GPI polarimetry (Monnier et al. 2017) is also not visible, likely due to heavy self-subtraction due to poor field rotation. The position of the Guidi et al. candidate lies well separated from the ring and residual speckle noise; the S/N maps show no convolved pixel within one PSF footprint (∼008) of this position with a significance greater than 1.3σ. More conservative reductions (e.g., larger rotation gap; higher SVD cutoff) may show slightly elevated residual emission consistent with additional extended structure at this separation (e.g., additional ring material). However, this signal is not statistically significant and is simpler to explain as residual speckle noise instead.

Figure 1 reveals the clear detection of scattered light surrounding the HD 163296 disk in our H-band imagery outside of the inner working angle of these data, 03 (30.5 au). The scattered light imagery reveals a broken ring structure that peaks at a distance along the major axis of 065 (66 au) and extends out to 098 (100 au) along the major axis (see Figure 4). Both the Q_ϕ and SN imagery exhibit little coherent signal between our inner working angle and the inner edge of the ring structure. We do not detect the inner disk component as previously detected by Monnier et al. (2017) due to our larger inner working angle. We therefore conclude that the small amount of scattered light interior to the ring in the PI image (panel A; Figure 1) could arise from a mixture of residual, uncorrected flux from the PSF and scattered flux from an inner disk component that is within our masked region (see, e.g., Takami et al. 2018). The NE side of the disk is known to be the near side (Rosenfeld et al. 2013), and the IR scattered light disk exhibits evidence of strong forward scattering (Guidi et al. 2018). The broken ring structure we observe is missing polarized intensity originating from the far side of the disk (SW region, along the minor axis; see Figure 1).

Zoom In Zoom Out Reset image size

We fit an ellipse to the scattered light ring using a least-squares fitting code (Hammel & Sullivan-Molina 2019),⁴² assuming the ring is a perfect circle projected at inclination. Since a known bias of the code is to prefer a smaller ellipse by preferably fitting the inner points (Halir & Flusser 1998), we choose to fit the peaks of the ring to mitigate this effect. Due to the low signal along the SW minor axis and sporadic structure along the NE minor axis, we did not keep any vertical cuts between 70° < PA < 180° and between 270° < PA < 370°. We fit a Gaussian to each vertical crosscut, producing the peak x, y position of the ring, and input these positions into the ellipse code described above.

In order to estimate the error of our ellipse fit, we performed a Monte Carlo routine by randomly sampling the Gaussian xy-coordinate errors and adding them to the xy-coordinates found above. We additionally applied a random rotation of the image between 0° and 1° to constrain the error associated with the interpolation of the image due to rotation. We performed 500 iterations and used the average values of the 500 iterations as the best-fit ellipse. The errors were estimated by taking the standard deviation of the parameters found with the 500 iterations. The best-fit results are shown in Table 1. The best-fit ellipse is compared with the PI image in Figure 5, shown as the white oval, along with the center of the disk (small white circle) and the center of the star.

Zoom In Zoom Out Reset image size

Our measured inclination of the disk (414 ± 03) and PA (1322± 03) is in agreement with the values derived from ALMA data of 42° and 132° respectively (Isella et al. 2016). Additionally, the offset of the minor axis from the central star that we find (−00432 ± 00016) is consistent with previous measurements, given their quoted errors when available (006, Garufi et al. 2014; 0105 ± 0045, Muro-Arena et al. 2018; 01, Monnier et al. 2017).

We applied an r² illumination correction to our data to better investigate the physical distribution of dust in the ring seen in Figure 1. We then azimuthally binned the average flux per area of the ring between two concentric ellipses. We adopted an inclination of 42°, and constructed each bin to be 8° wide and spanned a projected radial distance of 055–071 (55–72.5 au) to encompass the majority of the disk flux. The binned disk flux is azimuthally symmetric along the major axis, with the NW and SE side of the disks exhibiting the same amount of polarized intensity (Figure 6). There is also a clear azimuthal asymmetry in the binned flux along the minor axis, with the near side of the disk (NE side) exhibiting substantially more flux than the far side (SW side). We observe a deficit in scattered light flux along the near side of the disk at a PA of 30° in both the binned imagery (Figure 6) and unbinned PI, Q_rot images, and SN map images (Figure 1). This feature coincides with the position angle of the disk brightness enhancement and the position angle of the candidate point source noted by Guidi et al. (2018), and will be further discussed in Section 5.4.

Zoom In Zoom Out Reset image size

To help interpret our imagery and contemporaneous IR spectroscopy of HD 163296, we utilized the 3D MCRT code, HOCHUNK3D (Whitney et al. 2013). HOCHUNK3D allows the user to characterize the radial dust distribution, dust composition, and disk illumination parameters, and outputs a spectral energy distribution (SED) of the disk and imagery in a variety of user-defined bandpasses. The current version of HOCHUNK3D allows the user to decouple the disk into two dust distributions, allowing one to parameterize both a settled dust population toward the midplane and a different dust population in the upper surface layers of the disk. These two dust populations can either be co-spatial, or have different radial sizes. The dust distribution of each disk is characterized by several power-law parameters: the radial power law (α), the vertical Gaussian distribution (β), and the height of the disk from the midplane (h). Deviations from these power laws, such as a gap, spiral arms, warped disks, and walls, can all be included. The code also allows for the presence of a dusty envelope, which is parameterized by its minimum and maximum radius (${R}_{\mathrm{minenv}}$, ${R}_{\mathrm{maxenv}}$), and a dust density power law (ENVEXP). The dusty envelope can also include gaps and a bipolar cavity. Following the techniques established by Sitko et al. (2008), Wagner et al. (2015), and Fernandes et al. (2018), we use the dusty envelope as a proxy to model material ejected from the disk, aka a disk wind.

We constrained our model-starting parameters by observations when possible, and adopted the parameters from M. Pikhartova et al. (2019, in preparation), who are using HOCHUNK3D to model the variations seen in two epochs of HD163296's SED, as a starting point for our parameter-space exploration. ALMA observations of HD 163296 revealed the presence of 3 gaps located at 049, 082, and 131 (50, 83, and 133 au respectively given Gaia-DR2 distance of 101.5 pc) (Isella et al. 2016). Because our HiCIAO imagery is only sensitive to the first dust ring and the near-IR SED is most sensitive to dust features closer to the star, we only include the inner gap in our model. We allowed the two components of the dust distribution to be vertically stratified, and chose the radial extent of these distributions to match those observed for grains populating the midplane (250 au from VLA and ALMA observations; Guidi et al. 2016) and surface layers (540 au, Isella et al. 2007; Wisniewski et al. 2008). We note that while ALMA observations of the system were best described by a radial power law multiplied by an exponential function (Isella et al. 2016), HOCHUNK3D only uses a power-law function. Nevertheless, we did adjust the large grain dust distribution to match, as closely as possible, the dust distribution as measured by ALMA in the inner portion of the disk (Isella et al. 2016). The dust parameters for the large grain disk that we used are adopted from Wood et al. (2002), and are composed of amorphous carbon and silicon dust particles ranging in size up to 1 millimeter. The small grain disk and envelope dust parameters are from Kim et al. (1994), which is the average galactic interstellar medium (ISM) dust grain model.

We adopted an interstellar extinction of A_V = 0 mag from Ellerbroek et al. (2014), who measured the level of extinction from the ejected Herbig-Haro (HH) knots. Note that Ellerbroek et al. (2014) concluded that the optical variability of the SED likely comes from on-source reddening. In our model, we utilize the dusty envelope, a proxy to model disk wind, to replicate the on-source reddening, which is further discussed in Section 5.2. We explored accretion rates ranging from 1.73 × 10⁻⁷ to 4.35 × 10⁻⁶ M_⊙, calculated from contemporaneous Brγ emission lines in the SpeX 2011 data, first presented in Ellerbroek et al. (2014), but adjusted for the new distance of 101.5 pc.

We constrained these models using a SED (Figure 9) constructed from contemporaneous near-IR observations (Figure 2), along with non-contemporaneous photometry from the AllWISE catalog (Wright et al. 2010), 2MASS All Sky Survey (Cutri et al. 2003), IRAS point-source catalog (Helou & Walker 1988), and the historical variability of the V-band photometry as compiled in Ellerbroek et al. (2014). We also constrained these models using the surface brightness profiles along the major axis of our HiCIAO H-band scattered imagery (Figure 4).

We explored the parameter space of our models using a χ² minimization scheme. Namely, we calculated the χ² for the SED fit, the surface brightness along the major axis, and the minor axis offset, and added these values in quadrature to find the total χ² value. Because some of the SED data were not contemporaneous, we also calculated a separate χ² value that only incorporated comparisons of contemporaneously obtained data to the model. We began the iterative process with model runs of 5 million photons in order to find the best-fit SED to the SpeX and BASS spectra. Next, we increased the number of photons in each run to 50 million photons to obtain higher-quality model H-band images, and convolved the model image with the PSF of the H-band image. We explored parameter space to produce the best-fit χ² value between model surface brightness along the major axis and the minor axis offset to the observed imagery. After finding the best chi-squared fit model image, we iteratively switched between the SED and the model until we found a model that optimized the combined chi-squared value, resulting in our best-fit model. We then reran this best-fit model using 10⁹ photons to produce the model SED and imagery used in all of our figures. We remind readers that MCRT models, like HOCHUNK3D, that employ a large family of parameters suffer from parameter degeneracy; thus, our best-fit model is not unique (Dong et al. 2012).

Table 2 lists the main parameters utilized in our best-fit model, and Figure 7 details the temperature and density profile of the disk in this model. Figures 8 and 9 show the SED and radial surface brightness profile along the disk major axis of our best-fit model as compared with our observations. We remark that our best-fit model parameters are generally similar to those previously reported in the literature. For example, our disk mass of 0.05 M_⊙ (Table 2) is similar to that measured by Qi et al. (2011) (0.089 M_⊙) and Isella et al. (2007) (0.12 M_⊙).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  List of Key Best-fit Model Parameters and Estimates of the Upper and Lower Bounds of the Parameter

Parameter (Units)	Best-fit Model	Lower Bound	Upper Bound
Star Temperature (K)	9250	⋯	⋯
Star Radius (R⊙)	1.4	1.2	1.6
Disk Mass (M⊙)a	0.05	⋯	⋯
Fraction of Mass in Large Grain Disk	0.9	0.8	0.95
Inner Gap Radius (au)	29	20	32
Outer Gap Radius (au)	59	55	62
Large Grain Disk Minimum Radius (Rsub)b	31.9	25	35
Large Grain Disk Maximum Radius (au)	250.1	⋯	⋯
Large Grain Disk Scale Height (Rsub)b	0.11	0.08	0.13
Large Grain Disk Radial Density Exponent	0.1	0.05	0.2
Large Grain Disk Scale Height Exponent	0.16		0.18
Small Grain Disk Minimum Radius (Rsub)b	1.22	1.0	1.5
Small Grain Disk Maximum Radius (au)	540.1	⋯	⋯
Small Grain Disk Scale Height (Rsub)b	0.11	0.08	0.13
Small Grain Disk Radial Density Exponent	0.05
Small Grain Disk Scale Height Exponent	1.25
Envelope Inner Radius (Rsub)b	0.41	⋯	⋯
Envelope Outer Radius (au)	2.38	⋯	⋯
Envelope Density (${\rm{g}}\,{\mathrm{cm}}^{-3}$)	4.0 × 10−17	2.0 × 10−17	6.0 × 10−17
Accretion (M⊙)	6.0 × 10−7

Notes.

^aDisk mass value includes dust and gas. We assumed the gas to dust ratio is 100. ^bR_sub is the sublimation radius with 1 R_sub = 0.36 au.

Download table as:  ASCIITypeset image

Our best-fit model SED generally matches well with the contemporaneous spectroscopy and historical observations from optical to radio wavelengths (Figure 9). Since the optical flux has been shown to be highly variable and we do not have contemporaneous optical photometry or spectroscopy, we do not know whether the modest model overestimation of the optical flux simply reflects that the star was at a high flux state in 2011. Additionally, our model reproduces the on-source extinction value of A_V = 0.5 mag from Ellerbroek et al. (2014) with a value of A_V = 0.46 mag. We note that the observed versus model imagery comparison matches well along the NW side of the disk (right-hand side of Figure 8), while the model imagery is marginally too narrow along the SE side of the disk (left-hand side of Figure 8). This could be due to slight geometrical variations in the wall of the disk, causing the illumination of the SE side of the disk to be broader. We provide a full comparison of the observed H-band PI imagery and model imagery in Figure 10. Our model imagery reveals little scattered light beyond the bright ring and little to no scattered light within the gap of the disk, which matches the observed PI and Q_ϕ images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although none of the protoplanets/candidates reported from Keck/NIRC2 or ALMA are visible in our data, great care is needed to properly interpret these non-detections and their implications. For example, like HD 163296, HD 100546 has multiple imaged protoplanet candidates embedded in a bright, structured protoplanetary disk (Quanz et al. 2013; Currie et al. 2015). Follow-up claims of a spurious detection/non-detection of candidates around HD 100546 were faulty, as shown in Currie et al. (2017), in large part due to (1) incorrect spectrophotometric calibration and (2) a lack of forward-modeling of planet and disk signals.

Contemporaneous near-IR spectra of HD 163296 allowed us to spectrophotometrically calibrate CHARIS data cubes (see Section 2.1). To properly understand our non-detections and derive upper limits at the candidates' locations, we then performed forward-modeling of our images, investigating the reduction of the total source signal and the biasing of its spatial intensity distribution due to processing. This annealing results from self-subtraction of the source by itself and over-subtraction of the disk in ADI and SDI. Our method follows that outlined in Currie et al. (2018a), where we save the A-LOCI coefficients α and model the disk and planet signals as introducing a linear perturbation of value β, which provides an additional source of annealing (see also Brandt et al. 2013; Pueyo 2016). We focus on the 2018 May data due to its higher quality.

First, we explored the effect of the disk on the non-detections of planetary companions, using forward-modeling to determine its annealing due to processing and its effect on any point sources located exterior, like the proposed companions from Guidi et al. (2018) and Teague et al. (2018). We started with the best-fit scattered light disk model described in Section 3.2, which is drawn from our H-band scattered light imagery with Subaru/HiCIAO. We produced total intensity (not scattered light) images in J, H, and K passbands and interpolated the model images onto the CHARIS wavelength array and pixel scale. The model disk is slightly bluer than the combined light of the star+disk, with intrinsic colors of J–H, H–K of ∼0.35 and ∼0.35. Note that the model was constructed based on a single passband (H-band); thus, the model may not constrain the true color of the disk. The disk contrast with respect to the star on the forward-scattering side is typically ΔM = ${M}_{\mathrm{disk}/{\mathrm{arcsec}}^{2}}$ –M_⋆ ≈ 3.5–4. The visible trace of a disk may differ in total intensity versus scattered light. Therefore, we slightly adjusted the model parameters to provide a better match to the forward-modeled disk image, specifically increasing the semimajor axis by 5%.

Second, we verified that an object consistent with the 6–7 ${M}_{{\rm{J}}}$ candidate from Guidi et al. would be detected in our data. We used standard hot-start evolutionary models from Baraffe et al. (2003), adopting a planet age equal to the system age (5 Myr). This approach is intermediate between possible extremes that would yield higher and lower luminosities for a given planet mass. While we assume a planet age of 5 Myr when estimating mass limits, the age of a superjovian planet is likely much younger than that of the host star (Currie et al. 2013). This is especially true for protoplanets, which are nearing the end of their formation and thus much closer to t ≲ 1 Myr for any evolutionary model, where the planet luminosity is maximum. The inferred limits adopting would be then substantially lower than those we report. Conversely, we could adopt planet mass limits using the "cold-start" evolutionary models (e.g., Marley et al. 2007). However, recent literature casts serious doubt on the validity of the cold-start model formalism, which relies on specific assumptions about the entropy of accreted material. As shown by Berardo et al. (2017), classic cold-start conditions are extremely difficult to reach because the protoplanet will be substantially heated by the accretion shock, which will increase its entropy, resulting in hot-start-like initial conditions. Furthermore, imaged planets for which we have derived dynamical masses—β Pic b, HR 8799 bcde (Lagrange et al. 2010; Marois et al. 2010; Currie et al. 2011; Snellen & Brown 2018; Wang et al. 2018; Dupuy et al. 2019)—are inconsistent with a cold-start evolutionary model. At the candidate's location in each data cube, we injected a planet whose temperature matches that expected for a 4 M_J, 5 Myr planet according to these models. Although such a planet is predicted to be near the L/T dwarf transition (T_eff ∼ 1300 K), we assume a (cloudier) L dwarf spectrum drawn from the Bonnefoy et al. (2014) library, because annealing due to SDI will be stronger for such a spectrum. Integrated over the CHARIS wavelength array, the broadband contrast of this planet with respect to HD 163296 is ∼8 × 10⁻⁶, about 2.5–3.5 times as high as the predicted contrast for the Guidi et al. companion using a cloudy planet atmosphere from Currie et al. (2011).

Finally, our forward-modeling calculation allowed us to compute radially averaged, throughput-corrected broadband contrast curves. As with our fake planet injection test, we used the Baraffe et al. (2003) models to map between planet mass and temperature. To map between temperature and spectrum, we further used atmosphere models drawn from Burrows et al. (1997), adopting cloud prescriptions that provide reasonable fits to near-IR photometry for HR 8799 bcde and ROXs 42Bb, whose temperature, gravity, and masses cover most of our range (Currie et al. 2011, 2014; Madhusudhan et al. 2011; T. Currie et al. 2018, in preparation).

Figure 11 shows the wavelength-collapsed image of the input disk (left panel) and output image after forward-modeling the disk through ADI and SDI (right panel). While the disk in total intensity is more forward-scattering than the model based on polarimetry would predict, and its brightness is ∼30% higher, the model otherwise reproduces the CHARIS data and is sufficient for investigating the impact of self-subtraction on the forward-scattering side. The proposed candidate from Guidi et al. (2018) lies exterior to the main trace of the disk (cyan cross). After processing, the candidate's location is free of negative self-subtraction footprints. Inspection of the individual data cubes containing the disk model processed through ADI and SDI likewise show a flat background. At wider separations overlapping with the proposed candidate from Teague et al. (2018), the disk likewise leaves negligible residual effects.

Zoom In Zoom Out Reset image size

As shown in Figure 12 (left panel), a 6–7 M_J candidate similar to the one proposed in Guidi et al. should have been detected in our data. The fainter, even lower-mass (4 M_J) candidate injected into our data is clearly visible. While its S/N is formally ∼4.8, our inclusion of disk signal contributions leads our estimate of the noise to be conservative. A planet corresponding to the Guidi et al. candidate (ΔF ∼ 2.5 × 10⁻⁵) would be even more decisively detected (S/N ∼ 15).

Zoom In Zoom Out Reset image size

Broadband contrast limits in the right-hand panel of Figure 12 provide stringent limits on protoplanets covering the range probed with Keck/NIRC2 and ALMA. At ρ ∼ 049, the azimuthally averaged 5σ contrast limit is ∼8.5 × 10⁻⁶, in agreement with our expectations from the fake planet injection. If the Guidi et al. companion is real, it would then have to be redder than H–L_p ∼ 3.5 to escape detection: redder than all directly imaged planets except for the extreme L/T transition object HD 95086 b (De Rosa et al. 2016). Over the separations just interior or close to the visible trace of the disk and comparable to the separation of the Guidi et al. companion—ρ ∼ 04 (07) along the minor (major) axis—we can exclude planets with masses of 2–5 M_J, assuming standard hot-start evolutionary models. The CHARIS field encloses the possible location of the innermost companion proposed by Teague et al. (2018), which would lie at a projected separation of r_proj ∼ 83 au (ρ ∼ 082) along the major axis, or r_proj ∼ 40 au (ρ ∼ 04) along the minor axis. At these locations, our data rule out planets more massive than 5 M_J and ∼1.5 M_J, respectively. If located along the minor axis, the outermost proposed companion from Teague et al. (2018) would be at ρ ∼ 065 with a mass less than ∼2 M_J, according to our data.

HD 163296 has been observed numerous times across optical–IR bandpasses. Here we briefly summarize some of the major results of those investigations, to compare and contrast with our new imagery.

Space-based optical imagery has been obtained in both white light (HST/STIS; Grady et al. 2000) and broadband filters (HST/ACS; Wisniewski et al. 2008), tracing the disk out to 44 (447 au) and detecting HH knots. Comparison of these data revealed evidence for significant variation (∼1 mag) in the disk surface brightness, changes in the number of disk ansae visible over time, and changes in the relative brightness of features located in the NW and SE disk regions (Wisniewski et al. 2008). Unfortunately, none of these optical observations fully overlapped in wavelength coverage.

Ground-based AO imagery of the system can be generally summarized into 3 categories. First, a subset of observations clearly reveal the detection of the disk in scattered light, but the presence of residual AO speckle noise in the disk vicinity prevents a robust characterization of the surface brightness or detailed morphological structure of the disks (e.g., 2012 H-band imagery Garufi et al. 2014; 2014 Ks-band imagery Garufi et al. 2017). Second, a subset of observations (e.g., 2012 Ks-band imagery; Garufi et al. 2014) reveal the detection of an inclined ring structure extending out to 103 (103 au), where the intensity of scattered light is strongest along the major axis of the disk and is symmetrical about both sides of the disk major axis (NW side and SE side). Third, a subset of observations (e.g., 2014 J-band, Monnier et al. 2017; 2016 H-band imagery, Muro-Arena et al. 2018) reveal clear evidence of this same inclined ring structure whose flux is both azimuthally asymmetric and not the strongest along the major axis. In particular, the NW side of the major axis is brighter than the SE side of the disk in J-band GPI observations (see Figure 2, Monnier et al. 2017), and the maximum flux from the disk is north of the major axis peaking on the NW side of the disk in these data. The 2016 H-band VLT/SPHERE observations (Muro-Arena et al. 2018) also exhibit strong azimuthal asymmetry, with the NW side of the disk along the major axis exhibiting 2.7× more scattered light than the SE side of the disk along the major axis. Muro-Arena et al. (2018) used 3D radiative transfer modeling to suggest that this strong azimuthal asymmetry could be reproduced by including an inner disk component that was misaligned by 1° compared with the outer disk.

Our 2011 epoch H-band imagery is consistent with the second category of disk appearance we discussed in Section 5.1. Namely, we observe a broken ring structure in H-band scattered light whose flux peaks along the major axis and exhibits clear symmetry between the NW and SE side of the disk. Our 2011 epoch H-band data are thus clearly different from the 2016 epoch VLT/SPHERE H-band data, which show a 2.7× asymmetry between the NW and SE side of the disk (Muro-Arena et al. 2018).

To illustrate these differences, we scaled the peak flux along the major axis of the 2016 VLT/SPHERE data and present these data as dashed horizontal lines in our Figure 4. The 2.7× asymmetry about the major axis observed in the 2016 VLT/SPHERE data is clearly outside of the 3σ errors of our 2011 data. This obvious difference is also seen by comparing Figure 1 of Muro-Arena et al. (2018) with Figure 1 in this paper. We note that neither data set exhibits evidence of large-scale gradients in its U_ϕ component, indicating that systematic artifacts are not the cause of this phenomenon. We suggest that this is clear evidence that the system exhibits large changes in the appearance of its scattered light disk as seen in multi-epoch observations obtained with the same filter, and supports previous suggestions of this phenomenon as deduced from multi-epoch observations from similar, albeit not the same, filters (Wisniewski et al. 2008).

There are several mechanisms that could cause an azimuthal asymmetry of scattered light, including an asymmetrical distribution of dust (Muro-Arena et al. 2018), an inclined inner disk shadowing the outer disk (Muro-Arena et al. 2018), a warped inner disk structure shadowing the outer disk (Sitko et al. 2008), or dust ejected above the midplane of the disk that shadows the outer disk (Ellerbroek et al. 2014).

Muro-Arena et al. (2018) suggested that an asymmetric distribution of dust in the system was unlikely, as no asymmetry was observed with ALMA (Isella et al. 2016). Muro-Arena et al. (2018) were able to replicate the azimuthal asymmetry they observed in their scattered light imagery by inclining the inner disk by 1° compared with the outer disk, which is consistent with previous near-IR interferometric observations (Tannirkulam et al. 2008; Lazareff et al. 2017; Setterholm et al. 2018). However, our 2011 epoch data reveal the presence of no azimuthal asymmetry along the major axis in the same filter bandpass as the 2016 SPHERE observations. An inclined inner disk is unlikely to precess significantly over a 5 yr time frame; hence, an inclined inner disk alone is unlikely to produce the observed significant azimuthal variations in the scattered light disk. Moreover, we have shown that we can reproduce the basic properties of both our contemporaneously obtained near-IR SED and H-band imagery with a model that does not include an inner inclined disk. Thus, while the system could plausibly host an inclined disk, we suggest that this feature is unlikely to be responsible for producing the time-dependent azimuthal variations in the outer scattered light disk of the system.

We consider several other mechanisms that could explain the change in disk surface brightness seen in the system. First, a warped inner disk structure, such as a puffed-up inner disk wall (Turner et al. 2014), could be shadowing the outer disk (Sitko et al. 2008). If this disk warp were to dissipate or rotate azimuthally within a 2–3 yr timescale, this could cause a change in illumination of the outer disk similar to that observed between the 2011 and 2016 epoch H-band data sets. Dynamical simulations are needed to determine whether a substantial change in the appearance of a warped disk could occur on this short of a timescale and lead to the amplitude of variable disk illumination observed.

Second, this phenomenon could be caused by dust ejected above the midplane of the disk, which partially shadows the outer disk, as proposed by Ellerbroek et al. (2014). These dust "clouds" could differentially obstruct the illumination of the outer disk while they are between the star and the outer disk, as shown in Figure 13. We do have IR spectra that were obtained at a similar epoch to both our 2011 HICIAO data and the 2016 SPHERE data. The contemporaneous IR spectra cannot constrain the possible asymmetric nature of the dust clouds, but can constrain the total amount of dust in the disk wind when compared with our MCRT models. As shown in Figure 2, while both have the same flux around 0.9 μm, the 2011 epoch IR spectrum is brighter (∼0.5 mag at K') around 2 μm than the 2016 epoch IR spectrum. We remark that we can best reproduce the 2016 SED in our model by adopting a ∼2× lower envelope density, e.g., $9.0\times {10}^{-18}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, which corresponds to a lower circumstellar extinction in 2016 of A_V = 0.1 mag. We predict that the 2016 epoch should be 0.4 mag brighter in the V-band compared with the 2011 epoch data, similar to the optical variability found by Ellerbroek et al. (2014). Because the interstellar extinction is consistent with AV = 0 mag, the observed reddening most likely originates from the system. Thus an asymmetric disk wind launching dust clouds can explain both the variable illumination of the outer disk (Figure 13) and the reddening optical variability observed by our disk wind models and Ellerbroek et al. (2014).

Zoom In Zoom Out Reset image size

If the system does have an inclined inner disk as suggested by Muro-Arena et al. (2018) that during some epochs produces non-axisymmetric illumination of the outer disk (e.g., NW side brighter than SE side; 1998 HST/STIS, Grady et al. 2000, 2014 J-band, Monnier et al. 2017; 2016 H-band, Muro-Arena et al. 2018), the spatial distribution of any dust clouds elevated by a disk wind must also be non-axisymmetric to produce the observed epochs of axi-symmetric illumination of the outer disk (e.g., as seen in 2012 Ks-band imagery, Garufi et al. 2014; 2011 H-band, this study) and the sole-epoch of observed non-axisymmetric illumination with the SE side of the disk brighter than the NW side (2004 HST/ACS, Wisniewski et al. 2008). Future observations that simultaneously observe quiescent and wind events with contemporaneous optical and IR photometry and coronagraphic imagery could help to test whether shadowing by dust clouds could explain the observed behavior of the inner and outer disk of the system, and better parameterize the azimuthal distribution of such dust clouds.

We were able to reproduce the basic properties of our contemporaneous near-IR spectra and scattered light H-band imaging with a 3D MCRT disk model, which approximated the features of a disk wind via an envelope. As seen in Figure 9, our model SED is consistent with the highest observed V-band flux that was reported by Ellerbroek et al. (2014), but we caution that the robustness of this agreement is uncertain because we do not have contemporaneous optical photometry.

Muro-Arena et al. (2018) also performed MCRT modeling of HD 163296, and compared their models to the ALMA dust continuum image from Isella et al. (2016), their own VLT/SPHERE image, and historical photometry and spectroscopy. They modeled all three gaps that were observed in the ALMA continuum image and introduced an inclined disk to explain the asymmetric scattered light flux observed with the VLT/SPHERE image as noted above in Section 5.2. Their model images and SED are well matched to their observed images and historical photometry and spectroscopy. While they do not employ a disk wind model as we did (Section 3.2, see Figure 13), their model does not have a clear mechanism to explain the time-dependent azimuthal asymmetries seen in near-IR scattered light images (Section 5.2) or the optical–IR photometric and spectroscopic variability that has been observed (Sitko et al. 2008; Ellerbroek et al. 2014). We caution that the inability of an inclined disk by itself to explain the observed time-dependent azimuthal asymmetries in scattered light does not exclude the possibility that the system does in fact have an inclined disk. Due to limitations with HOCHUNK3D, we leave applying our disk wind model to the archival images and SEDs to future work.

We note that a deficit of scattered light is seen in the near side of our disk imagery at a PA of 30°, in both our binned imagery (Figure 6) and our unbinned PI and Q_rot images. We caution that while this feature could be real, it is not uncommon to observe depolarization along the minor axis due to the residual presence of an uncorrected polarized halo. Interestingly, this feature coincides with the disk brightness increase observed with the Keck/NIRC2 L'-band vortex coronagraph by Guidi et al. (2018) and is located at a similar position angle, albeit closer to the host star, as the purported candidate planetary mass object reported by Guidi et al. (2018). As noted by Guidi et al. (2018), this disk feature is located where forward scattering should be significant. If the feature we observe at the similar disk position is astrophysical, the decreased amplitude of the feature in polarized intensity suggests that it could be polarized less than its neighboring disk material.

Our data improve the detection limits for protoplanets in thermal emission around HD 163296 compared with Keck/NIRC2 data from Guidi et al. (2018): from 5–7 M_J to now 2–5 M_J, assuming standard hot-start evolutionary models, near the projected trace of the disk. At wider separations covering the possible locations of the inner proposed candidate from Teague et al. (2018) (r_proj ∼ 83 au/ρ ∼ 082), the limits have now improved from 4.5 M_J to 1.5 M_J, the latter of which is just slightly higher than the predicted mass of the companion (1 M_J). Limits for the outer Teague et al. candidate along the minor axis are likewise just slightly higher than the predicted mass (a limit of 2 M_J versus a predicted 1.3 M_J). Thus, at least for now, the ALMA-predicted protoplanet candidates are consistent with direct imaging constraints.

Our data appear to rule out the proposed, marginally significant candidate identified from thermal IR data in Guidi et al. (2018). Using standard assumptions for planet atmospheres, our forward-modeling demonstrates we could have detected an even fainter planet at the location of the proposed candidate. For an assumed age of 5 Myr and hot-start evolutionary models, the candidate is predicted to be 6–7 M_J, while our radially averaged contrast limits are significantly lower (∼4–5 M_J).⁴³

The simplest explanation for our conflicting results is that the NIRC2 candidate is instead residual, partially subtracted speckle noise or partially subtracted disk emission left over from processing. Figure 1 of Guidi et al. (2018) shows multiple emission peaks with a similar or slightly smaller spatial scale as the candidate (e.g., at the 2, 6, 7, and 8 o'clock positions just exterior to the masked region). An even brighter, seemingly point-source-like peak at nearly the same position angle in these data appears to be an artificially enhanced region of the disk, which could have been mistaken for a point in shallower and/or higher background data. Convolving the image with a Gaussian kernel may further accentuate the point-source-like appearance of these features.⁴⁴ The position of the candidate also coincides with the minor axis of a second ring of emission detected with ALMA. Forward-modeling as performed in Currie et al. (2015) could better clarify whether the candidate's morphology is consistent with an annealed point source or residual disk emission.

Alternatively, the candidate could be extremely red/underluminous in the near-IR and thus difficult to detect. If embedded in the disk, it would be preferentially extincted in the near-IR compared with the thermal infrared, as has been proposed for HD 100546 b (Currie et al. 2015; Quanz et al. 2015). It could also retain an extremely dusty/cloudy atmosphere characteristic of some young exoplanets near the L/T transition (Currie et al. 2011; De Rosa et al. 2016), making it appear "underluminous" in the near-infrared. Follow-up thermal infrared imaging at L_p or M_p could provide a more decisive probe of these possibilities.