Complex Spiral Structure in the HD 100546 Transitional Disk as Revealed by GPI and MagAO

Initial GPI observations of HD 100546 were taken in H-band spectroscopic mode (hereafter H-spec) using Angular Differential Imaging (ADI; Marois et al. 2006) as part of the Gemini Planet Imager Exoplanet Survey (GPIES). Follow-up observations were conducted in both spectroscopic and polarimetric modes based on extended structures suggested by this preliminary data set. A full summary of the GPIES observations is given in Table 1. All initial reductions were done using the GPI Data Reduction Pipeline (DRP) version 1.3.0 (Perrin et al. 2014, 2016). We refer the reader to these papers for full details of the GPI DRP. In brief, the GPI DRP subtracts dark background, interpolates over bad pixels, corrects for DC offsets in the 32 readout channels, and converts from raw 2D IFS frames to 3D datacubes. In the case of spectral data, argon arc lamp exposures taken both at the beginning of the night and immediately prior to the science exposure sequence are used for wavelength calibration, and the locations and fluxes of the four satellite spots created by the apodizer are used to compute and apply astrometric and photometric calibrations. In the case of polarimetric data, the pipeline assembles a full Stokes datacube from the sequence of exposures.

Table 1.  Summary of Gemini and Magellan Data Sets

Instrument	Date	Observing Mode	nimages	tint/frame	ncoadds	tint total	θrot	Avg. seeinga
				(sec)		(min)	(deg)	('')
MagAO	2014 Apr 11	Hα SDI	3423	2.273	1	129.7	73.5b	1.05
MagAO	2014 Apr 12	Hα SDI	4939	2.273	1	187.1	71.6	0.58
GPI	2014 Dec 17a	H-spec	33	60	1	33	12.9	⋯
GPI	2015 Jan 30a	Y-pol	14	60	1	14	⋯	0.63
MagAO	2015 May 15	Hα SDI	2077	2.273	1	78.7	42.0	0.46
GPI	2016 Feb 27	H-spec	120	60	1	120	51.6	0.66
GPI	2016 Jan 28	Y-pol	62	15	4	62	⋯	0.69

Notes. Bolded rows represent the three highest-quality data sets, which are used for the bulk of the analyses in this paper.

^aThe instrument and method for measuring seeing varies by telescope and observing run. Magellan seeing values are derived from measurements taken at the Baade telescope. Gemini South has both MASS and DIMM seeing monitors; however, only the DIMM was functioning on 2014 April 11, 2014 April 12, and 2015 January 30 and neither was functioning on 2014 December 17. On 2015 February 27 and 2016 February 28, both were online and seeing recorded by the two instruments has been averaged. ^bAlthough this data set has a total of $73\buildrel{\circ}\over{.} 5$ rotation, the space is not evenly sampled and there is a 10° gap in rotational space while the system was pointed at an NIR reference PSF star.

Download table as:  ASCIITypeset image

2.1.1. Polarimetric Data

Y-band polarimetric images (hereafter Y-pol) were first attempted on 2015 January 30; however, the sequence was aborted due to poor conditions. The data set we analyze in this work was collected on 2016 January 28. Data were taken using the shortest GPI filter (Y-band, 0.95–1.14 μm) and accompanying Y-band coronagraph because this mode affords the highest angular resolution and has the smallest coronagraphic inner working angle (0078), allowing us to probe the very inner regions of the disk near the HD 100546 "c" source and inside the inner cavity rim. Standard PDI waveplate cycles of $0^\circ ,22\buildrel{\circ}\over{.} 5,45^\circ ,$ and $67\buildrel{\circ}\over{.} 5$ were employed in order to allow for double difference polarized imaging (Kuhn et al. 2001; Perrin et al. 2004; Hinkley et al. 2009; Hashimoto et al. 2011), and the 60 s images were taken as four 15 s coadds in order to avoid saturating the inner region of the images.

Generation of the Stokes cubes was done using the standard polarimetry recipes available in the GPI pipeline and described in detail in Perrin et al. (2014) and Millar-Blanchaer et al. (2015); however, three modifications were made to the standard polarization recipe. First, we assumed a perfect half waveplate rather than the default lab-measured waveplate retardance. We have found in some cases that this improves the final image quality, as measured by the amount of residual signal in the U_r image (see next paragraph). The mean stellar polarization was estimated for each datacube individually from the normalized difference of the two orthogonal polarization slices in the region $0\lt r\lt 5$ pixels ($r\lt 0\buildrel{\prime\prime}\over{.} 07$), which is beneath the focal plane mask. Light that lies in this region should be composed primarily of starlight diffracted around the FPM, and any polarized signal is most likely induced by the instrument optics if we assume that the starlight is intrinsically unpolarized. The mean normalized difference in this region is scaled to the total polarized flux in each pixel before removal. For more details about the specifics of this estimation, see Millar-Blanchaer et al. (2016). Finally, we smoothed the processed images with a two-pixel FWHM Gaussian kernel before combining into the Stokes cube (I, Q, U, V) in order to mitigate microphonics noise (Ingraham et al. 2014).

The Stokes cube generated by the GPI DRP was transformed to a radial Stokes cube (I, ${Q}_{\phi },{U}_{\phi }$, V; see Schmid et al. 2006) via the same method as in Millar-Blanchaer et al. (2015). Under this convention, all polarized signal oriented parallel or perpendicular to the vector connecting the pixel to the central star is encompassed in the ${Q}_{\phi }$ image, and all signal oriented at ±45° is encompassed in the ${U}_{\phi }$ signal. The ${Q}_{\phi }$ image thus contains the centrosymmetric polarized disk signal in the case of single scattering, and the ${U}_{\phi }$ image is an approximation of the noise, under the assumption that the contribution of multiple-scattered photons is small.

The final ${Q}_{\phi }$ and ${U}_{\phi }$ images are shown in Figure 1 and discussed in detail in Section 3.1.

Zoom In Zoom Out Reset image size

2.1.2. Spectroscopic Data

The H-spec coronagraphic data set taken on 2014 December 27, while nominally a full GPI sequence, had low overall field rotation ($12\buildrel{\circ}\over{.} 9$), and its utility was compromised as a result. For this moderate inclination and highly extended disk, a large amount of field rotation is necessary to minimize disk self-subtraction and extract robust disk structure, therefore we followed up this initial observation on 2016 February 27 with a two-hour on-sky sequence, reaching $51\buildrel{\circ}\over{.} 6$ field rotation under good weather conditions (see Table 1). Results based on these later observations are presented in this paper.

The GPI DRP processes the raw data through dark subtraction, wavelength calibration based on observations of an argon arc lamp (Wolff et al. 2014), bad pixel identification and interpolation, microspectra extraction to create $(x,y,\lambda )$ datacubes (Maire et al. 2014), interpolation to a common wavelength axis, and distortion correction (Konopacky et al. 2014). Astrometric calibration (platescale of 14.166 ± 0.007 mas/pixel, position angle offset of $-0.10\pm 0\buildrel{\circ}\over{.} 13$) was obtained with observations of the ${\theta }_{1}$ Ori field and other calibration binaries following the procedure described in Konopacky et al. (2014).

Further post-processing was also done using the GPI DRP. The 3D datacubes were first aligned using the photocenter of the four satellite spot positions (Wang et al. 2014). To remove slowly evolving large-scale structures, the datacubes were high-pass filtered using a smooth Fourier filter with cutoff frequencies between 4 and 16 equivalent-pixels in the image framework, allowing us to investigate disk features on different spatial scales. Since this step strongly affects the apparent geometry of the disk, the two extremes of these cutoff frequencies, as well as images without any high-pass filter applied, are discussed and shown in Section 3.2.

The stellar PSF was estimated and subtracted from each image in the sequence using several ADI algorithms: classical Angular Differential Imaging (cADI; Marois et al. 2006), Locally Optimized Combinations of Images (LOCI; Lafrenière et al. 2007), and Karhunen–Loève Image Processing (KLIP, a form of Principal Component Analysis; Amara & Quanz 2012; Soummer et al. 2012) via a custom IDL pipeline. Using different ADI algorithms was mandatory in the analysis of this inclined, asymmetric, bright, and extended transitional disk to better assess the robustness of resolved structures against residual speckles, which manifest themselves differently in the post-processed images computed by each algorithm. For all three algorithms, residual images were rotated to align north with the vertical, combined with a 10% trimmed mean (discarding the highest and lowest 5% of pixel values in the temporal sequence), and collapsed over the wavelength axis to create a final broadband image.

cADI processing has no tunable parameters. The stellar PSF subtracted from each image is simply the median of the entire image cube, and the PSF-subtracted images are then rotated to a common on-sky orientation before combining. This method therefore is not capable of removing evolving PSF features, but it provides a good estimate of the most static PSF structures. Though mitigated by the large amount of field rotation, the disk extends azimuthally over more than the 516 of rotation in the data set, so some disk emission survives into the median PSF, resulting in negative "self-subtraction" regions at the edges of the disk.

LOCI analysis was done with annuli of dr = 5 pixels, optimization region of N_A = 500 FWHM (3.6 pixels at the H-band), geometry factor g = 1, and minimum separation criterion ${N}_{\delta }=1$ FWHM.

KLIP analysis was done on a single image region from 5 to 100 pixels in radius (007–142) and keeping only the first one to five Karhunen–Loève (KL) modes. Although KLIP is typically used for point-source searches with more zones and a greater number of KL modes, this single-zone, small number of KL-mode approach is standard for minimizing self-subtraction of extended disk features.

The PSF was also subtracted using the Reference Differential Imaging technique (RDI) implemented in the TLOCI quick-look processing pipeline (an evolution of the SOSIE pipeline; Marois et al. 2010). A library of reference images was created from 426 H-band datacubes (all GPIES campaign observations taken in pupil-stabilized mode at the H-band at the time of processing). Data from each reference sequence were first reduced with the GPI DRP in the standard manner described previously. Additionally, each image in an object sequence was high-pass filtered using an 11 pixel (016, 4 λ/D) square unsharp mask, magnified to align speckles across wavelength channels, flux normalized so that the satellite spot intensities were the same in each channel, and wavelength collapsed (only slices 5–31 were used to avoid the noisy wavelength slices at the end of every GPI spectral cube). These high-pass filtered, aligned, normalized, and wavelength-collapsed images were then median-combined for each object sequence and scaled to the flux of the target star using the satellite spots, allowing us to gather a homogeneous library of achromatic speckle-limited images with greatly reduced disk, planet, or background star signals. The HD 100546 DRP images were processed through the TLOCI RDI pipeline using only the 20 most correlated reference images in this PSF library to subtract the speckle noise. Reference images were selected by performing a cross-correlation analysis in a [15–80] pixel (5.4–28.6 $\lambda /D$, 0212–1133) annular region to avoid the focal plane mask edge.

RDI, cADI, LOCI, and KLIP images are shown in Figure 2 and are discussed in detail in Section 3.2.

Zoom In Zoom Out Reset image size

High-contrast, visible light, AO observations were conducted at the Magellan Clay telescope with the MagAO System (Morzinski et al. 2016, 2014) and its visible light camera (VisAO; Males et al. 2014). The observations were conducted in Hα Simultaneous Differential Imaging (SDI) mode in which a Wollaston prism is used to split the beam into two channels, and each is passed through a separate narrowband filter, one centered on the Hα emission line (656 nm, Δλ = 6 nm) and one centered on the nearby continuum (642 nm, Δλ = 6 nm). The continuum image serves as a sensitive and simultaneous probe of the stellar PSF.

GAPlanetS data are reduced with a custom IDL pipeline. Raw data frames are bias subtracted and divided by a flat field generated from R-band twilight sky observations. Dust spots in the instrument optics create significant throughput effects and can create point-like artifacts, but are clearly revealed in the flat field. They are not effectively removed by simply dividing by the flat field (why is unclear), and we therefore mask all pixels within 2 pixels of a region with <98% throughput. This mask is applied to all data frames before further analysis.

The bias-subtracted and flat-fielded raw images are then separated into line (e.g., Hα) and continuum channels. Individual channel images are registered against a high-quality individual data frame. The center of rotation is found through a custom centering algorithm that locates the center of circular symmetry in the median collapsed registered datacube by finding the point that minimizes the standard deviation of intensity in annuli centered at that point. We find that this algorithm performs better than radon transform or center of rotational symmetry algorithms for VisAO data, as measured using a binary of well-known separation and P.A.

Although a minimum integration time of 2.273 s was used in all cases, the HD 100546 observations were saturated at radii interior to 7 pixels in all data sets. This was noted during observations; however, we opted to saturate the very inner region rather than decrease detector gains and therefore sensitivity. We apply a software mask of r = 8 pixels to all data before centering, and exclude pixels interior to this radius in our KLIP reductions.

It is important to note that although Hα emission is typically thought of as an accretion tracer, disk-scattered light also makes a significant contribution at this wavelength, particularly in cases where the star itself is actively accreting. In fact, disk-scattered light is ubiquitous even in PSF-subtracted images because of the extended and moderately inclined nature of the HD 100546 disk. The MagAO system was designed to utilize the simultaneous nature of our Hα and continuum observations to remove both direct starlight and disk-scattered light contributions that are equivalent in the two filters. We compensate for the difference in stellar (and therefore scattered light) brightness between the two filters by subtracting a scaled version of the simultaneous continuum image from each Hα image before further processing. The scaling factor is determined iteratively as the value that results in minimized noise residuals in the region $8\lt r\lt 27$ pixels (representing the region between where the saturated images reach linearity and the inner boundary of the AO control radius). Scaling and subtracting the continuum image in this way effectively removes the contribution of scattered-light disk structures and diffracted starlight from the images. Because accreting protoplanets are expected to exhibit Hα excess and to not have a detectable level of continuum emission, this strategy should eliminate starlight and disk-scattered light preferentially, leaving behind pure Hα emission. It is the KLIP-processed versions of these SDI images that we use to place constraints on Hα emission from accreting protoplanets in these data sets.

KLIP images are generated using the MagAO interface of pyKLIP, a Python implementation of the KLIP algorithm (Wang et al. 2015). Of particular importance to the discussion in this paper is the fact that the final images are very sensitive to our choice of KLIP parameters, notably zone size and masking parameters. Although not exhaustive, we explore a wide region of this parameter space in order to assess the robustness of the parameters we extract, as reported in Section 3. The AO control radius for the MagAO system lies at r = 35 pixels in our images, and we find that the region $27\lt r\lt 42$ pixels is particularly noisy as a result, with many short-lived speckles that are not well-subtracted with KLIP. We therefore mask this region in each image before KLIP processing to avoid the appearance of spurious structures in the final reductions.

The moderate inclination of the disk means that disk structures cover wide swaths in azimuth, and aggressive KLIP reductions can be potentially problematic. We find that KLIP reductions with small to moderate exclusion criteria (e.g., allowing images where a planet would have moved by fewer than 8 pixels within a given annulus) result in large heavily self-subtracted regions and turn extended disk features into spurious point sources.

We elect the least-aggressive exclusion criterion possible for each data set, excluding all images from the reference library where a hypothetical planet located in the center of an annulus would have moved by fewer than a given number of pixels, where that number is as large as possible. For the 2014 data set, this is 12 pixels, corresponding to 33° of rotation in the innermost annulus before an image is included in the reference library. For the 2015 data set, this is 8 pixels, corresponding to 21° of rotation. Since the maximum exclusion criterion is nearly twice as aggressive in the case of the 2015 data set, it is unsurprising that the disk rim is not as cleanly revealed as in the 2014 reductions.

KLIP-processed Hα, Continuum, and SDI images for both epochs are shown in Figure 3 and discussed in detail in Section 3.3.

Zoom In Zoom Out Reset image size

GPI Y-pol images, shown in Figure 1, clearly resolve the scattered light cavity rim. There are distinct bright lobes along the disk major axis; however, these are symmetric about the star, and we see no evidence in these data of anything unusual at the location of the purported HD 100546 "c" point source.

The corresponding GPI ${U}_{\phi }$ image shows a non-zero signal, peaking at ∼20% of the value of the ${Q}_{\phi }$ image with most of the signal localized east of the star and just outside the coronagraph. This is potentially an effect of instrumental polarization, but a non-zero signal in ${U}_{\phi }$ HD 100546 images has been seen before (albeit with a different signal morphology; Avenhaus et al. 2014; Garufi et al. 2016), and may be a result of physical rather than instrumental effects. For example, multiple scattering is expected to create non-zero ${U}_{\phi }$ signals (Canovas et al. 2015).

In order to compensate for the purely geometric ${r}^{-2}$ dropoff in stellar scattered light, we scaled the images by r² for a disk inclined at 42° along a P.A. of 145°, a common practice in the field for revealing fainter extended structures in the outer disk. We note that we apply this scaling only to highlight faint disk features and that any asymmetries in brightness or location of disk features along the minor axis are impacted by the inclined, vertically extended and optically thick nature of the disk, which will tend to artificially enhance the illuminated half of the disk. The r²-scaled images do, however, effectively reveal a faint extended feature connected to the southeastern disk rim and extending to the north, which we will refer to hereafter as the "northeastern spiral." This feature is also effectively revealed with a simple four-pixel Fourier high-pass filtering of the original image. This and other morphological features revealed in GPI and MagAO imagery are discussed in detail in Section 4.1.

Radial profiles taken through the GPI Y-pol images, shown in Figure 4, reveal that there is no significant deviation between profiles taken to the east and west along the major axis, despite the proposed existence of a planet candidate along the eastern major axis. The profiles peak at 014, suggesting a cavity rim at 15 au. This is marginally inconsistent with the cavity radius estimated with SPHERE at an R of 12.5 ± 1 au, but quite consistent with the range of estimates (15–17 au) in the literature for the NIR cavity rim.

Zoom In Zoom Out Reset image size

The minor axis profiles are significantly different both in radial extent and in absolute intensity along the northern and southern minor axes; however, this is an expected effect. The greater radial extent and brightness of the northern minor axis profile is consistent with it being the illuminated half of the disk, and is likely affected by both the geometry of the disk and the scattering phase function.

Given the dearth of successful detections of polarized light from young planets in the literature (only upper limits; e.g., Jensen-Clem et al. 2016), it is perhaps unsurprising that there is no evidence of a point source at the location of HD 100546 "c" in the Y-pol image; however, there is a clear polarized disk structure at this location, and its smoothness and symmetry with respect to disk features opposite the star are surprising in the context of a planet at or near this location.

PSF-subtracted H-band images processed through a variety of reduction techniques are shown in Figure 2, and these techniques increase in aggressiveness toward the bottom of the figure. The apparent morphology is somewhat sensitive to the image processing technique. In particular, less aggressive PSF subtraction techniques (RDI, cADI) result in images that are dominated by an arc of emission extending from SE to NW. A number of additional, fainter structures resembling spiral arms are present to the south and east of the star, including several in the RDI and cADI images. Aggressive processing with LOCI and KLIP highlights these features further and reveals additional fainter structures; however, these aggressive techniques suppress the more extended arc of emission apparent in the cADI and RDI reductions.

MagAO images are shown for both the Hα and continuum channels, as well as SDI images (Hα−scale × continuum) in Figure 3. The structures in processed continuum images closely mimic the structures in the Hα images, which point to their common origin as disk-scattered light. Both images reveal an arc of emission consistent with the forward-scattering portion of the disk rim.

The importance of field rotation in the identification of high-fidelity disk features is apparent in the 15 May 2015 images, which had significantly less field rotation (420) than the 12 April 2014 images (716). The same forward-scattering inner disk rim is seen in this case, but it appears clumpy, and structures along it might even be mistaken for point sources.

The SDI images for both data sets, on the other hand, are free of extended scattered light structures. This points to the effectiveness of the process of scaling and subtracting the continuum image before KLIP processing. The images are also, unfortunately, free of any Hα excess point source candidates. This is perhaps unsurprising at the location of the b candidate, as it is embedded in the disk and very little dusty material is needed to extinct visible light emission. However, it is somewhat surprising at the location of the "c" planet candidate, which should be minimally extincted if it lies interior to the disk rim and inside of the relatively dust-free disk cavity. Quantitative constraints on detectable contrast levels for the "c" planet are discussed in Section 4.2.

It is important to note that the MagAO images presented here have markedly lower Strehl ratios than the GPI images, due to the fact that AO correction is significantly more difficult to accomplish in the visible than in the NIR since a given optical path difference will correspond to a larger fraction of a wavelength in visible light and naturally produce lower Strehl ratios. How much lower is difficult to estimate given the difficulty of measuring Strehl ratios in general and in saturated data in particular, but they are on the order of ∼10%–20% with MagAO at Hα, ∼25%–35% for the GPI Y-pol data set, and ∼65%–75% for the GPI H-spec data set. At the same time, the MagAO images benefit from the higher resolution afforded by visible light imaging, which compensates in part for the lower Strehl imagery.

In the previous section, we discussed features revealed in each data set individually. Here, we discuss how these multiwavelength data complement one another. With the exception of polarized data, where post-processing is minimal, it is unclear from a single data set alone whether all apparent disk features are true disk structures or artifacts of overly aggressive PSF subtraction processes. Such techniques have two problems when applied to disks in general, and moderately inclined disks for which significant disk emission survives into reference PSFs in particular. First, surface brightness measurements are severely complicated by disk self-subtraction (Milli et al. 2012), and we therefore do not attempt them in this work. Second, the morphology of complex disk structures can be compromised and spurious point-like artifacts introduced by self-subtraction. By overlaying the three data sets we have obtained and comparing them with features identified previously in the literature, we attempt to address this second point and identify the most robust disk features.

All three data sets can be seen on the same angular scale in the top panel of Figure 5, and the bottom panel shows pairs of images overlaid on one another. The smaller coronagraphic mask in the GPI Y-pol data set and the non-coronagraphic MagAO data allow us to fill in features in the very inner disk region, and the higher sensitivity of the H-spec data allows us to probe features in the outer disk. We have selected the H-spec cADI processed data set with a four-pixel high-pass filter for this analysis as it is less aggressive than the KLIP and LOCI-processed images, but reveals more of the faint disk features than the other cADI images (the high-pass filter serves to sharpen the disk features and therefore mitigates the azimuthal extent of the self-subtraction). The overlays reveal several very robust features present in multiple data sets, including the inner disk rim and the northeastern spiral arm.

Zoom In Zoom Out Reset image size

We label the most prominent revealed features from all three data sets in Figure 6 and discuss them below. We aim here simply to identify and name the most robust features and to compare them to features previously identified in the literature. A detailed discussion of the physical nature of these features, and the spiral arms in particular, is beyond the scope of this work, although we do engage in a brief qualitative comparison with spiral disk models viewed at moderate inclination in Section 4.3.

Zoom In Zoom Out Reset image size

Global Near/Far Side Asymmetry—The near side of the disk (inclined toward the observer, here the SW side) appears mostly featureless in all three images, whereas most of the structures are present on the far side (NE). This is a natural effect of observing an inclined flared disk, wherein the near-side disk geometry causes surface features to be compressed in projection or even shadowed from view by the disk midplane. The H-spec data also reveal a bright lane to the southwest, indicated with an "A" in Figure 6. This feature may be the front edge of the bottom (opposite the disk midplane) side of the disk, as discussed in Section 4.3. A similarly offset bright lane feature was recently detected by de Boer et al. (2016) in the disk of RXJ1615.3-3255 (Feature A1).

Inner Cavity—The inner cavity rim seen in both our GPI Y-pol data (Figure 1) and MagAO data (Figure 3) and indicated with a cyan ellipse in Figure 6 is extremely robust. Its existence is consistent with the NIR deficit in the SED of HD 100546 and with previous resolved images with VLT/NaCo (Avenhaus et al. 2014) and VLT/SPHERE-ZIMPOL (Garufi et al. 2016), though its location in the Y-pol radial profiles is marginally inconsistent with the latter. The potential for disk self-subtraction to affect the apparent location of the disk rim, as well as the close proximity to the H-band coronagraph preclude robust measurement of the disk rim location in total intensity at the H-band or Hα. Therefore, we defer discussion of whether the marginal inconsistency of our Y-pol disk rim radius with the shorter-wavelength SPHERE data is a wavelength-dependent effect for future work.

Disk "Wings"—All three of our data sets also reveal an extended arc of emission that runs through and beyond the southern rim of the disk cavity. With aggressive processing, this rim feature can appear sharp, but less aggressive subtractions suggest that it is in fact quite extended. It coincides with the sharp features labeled S5 and S1 in Figure 6, but can best be seen in its extended form in the cADI and RDI images of Figure 2. It is unclear whether the sharper features that we have labeled S5 and S1 are spirals embedded in that bright wing of emission or are that same feature made sharper by ADI processing. These "wing" features are the brightest and most distinct features far from the star, and have been identified in several previous studies (Currie et al. 2014, 2015; Garufi et al. 2016).

Spiral Arms—The spiral feature labeled S3 in Figure 6 is clearly visible in the minimally processed Y-pol data, and this also coincides with a brighter region in the Magellan data, though only a portion of it is visible inside of the masked AO control radius region, as revealed by the lower right panel in Figure 5. The Y-pol structure is also contiguous with H-spec emission that curves toward the feature labeled S2, and it is likely that these two features are part of the same spiral arm. This S3–S2 arm was also seen, though similarly broken, in the deep SPHERE/ZIMPOL polarimetric imagery reported in Garufi et al. (2016) and by Avenhaus et al. (2014). The geometry of this feature is puzzling if it is contiguous, as the apparent curvature back toward the star would suggest that S2 is at least in part a near-side feature, yet it does not obscure the cavity. Future deep polarimetric imaging is needed to understand the nature of this feature.

The inner parts of the S1 and S5 features are coincident with the disk "wings" described above, but the S1 feature curves inward more sharply and is consistent with the "Northern arm" identified in Garufi et al. (2016). It may be contiguous with the feature labeled S4, though, like the apparent S3–S2 spiral, this S4–S1 spiral is broken. The S4 spiral feature is faint and lies in a region near the bright disk wings that is especially heavily affected by disk self-subtraction, but it too has been seen in previous imagery and is labeled "spiral 2" in Currie et al. (2015).

The spiral feature S6 is also apparent in both MagAO and GPI H-spec data, though there is a break in the revealed feature approximately midway along the line labeled S6. This is the only such feature present on the near (SW) side of the disk major axis in our data. It may be a continuation of a spiral originating on the other side of the disk (S3/2 or S4/1), or it may be a secondary spiral arm mirroring a northern spiral. Similar "Southern Spirals" were identified in Garufi et al. (2016), albeit farther out. The Garufi et al. (2016) SPHERE K-band total intensity images reveal the same feature we have identified as S6, though it is not labeled by the authors as a feature of particular note.

We engaged in a brief exploratory modeling effort, described in Section 4.3, in an attempt to understand the identified disk structures and the effect that PSF processing can have on them. However, much work remains to be done in this area.

Our MagAO Hα and continuum data and GPI H-spec data reveal a bright apparent point source at r = 145 mas and P.A. = 152° after aggressive PSF subtraction and/or aggressive high-pass filtering (see Figure 7). This is consistent with the location of the candidate protoplanet put forward by Brittain et al. (2014) and supported by analyses in Currie et al. (2015) using a previous GPI H-band data set, so appears at first glance to be a promising planet candidate detection.

Zoom In Zoom Out Reset image size

However, it can be seen that this apparent point source is located at the intersection of the disk rim with the northeastern spiral arm and is mirrored by another concentrated knot of emission on the opposite side of the major axis. Although the symmetry of these features and coincidence with spiral arm intersections do not definitively rule out the existence of an underlying point source at the location of the "c" candidate, they do raise questions regarding its nature. The discovery paper by Currie et al. (2015) allowed for the possibility that this feature is a disk artifact, and we explore that scenario in this section.

To assess the hypothesis that the "c" candidate is a disk artifact, we engaged in two lines of inquiry.

GPI Spectra—The contrasts of two knots of emission (indicated with circles in Figure 7), one at the location of the "c" candidate and the other at the same location on the opposite side of the star, were extracted from our GPI H-band data using aperture photometry with a radius of 0.75 × FWHM (3.6 pixels) using the four-pixel high-pass filtered PCA (KL = 1) reduced wavelength images. Spectra of these knots were obtained after normalization with the spectrum of the star, obtained from the average of 10,400 and 10,600 K BT-NextGen models (Allard et al. 2012) and binned to the resolution of GPI. Since the two knots lie at the same stellocentric separation, they suffer from equivalent self-subtraction due to ADI and so have the same approximate uncertainties. Since we were only interested in the ratio of the two spectra, we rely on this symmetry to cancel out systematics due to PSF subtraction processing. Results are shown in Figure 7 (right panel). Not only does the spectrum of the source at the location of candidate "c" closely match the spectrum of the opposing knot of emission, it also shows no significant deviation from the spectrum of the star, pointing to a scattered light disk origin and showing no indication of an underlying planetary photosphere.

MagAO SDI Imagery—If the "c" candidate were indeed a protoplanet lying inside the disk gap, we might expect it to be actively accreting as gas passes through the dust cavity en route to the still-accreting central star. The cavity is also depleted in small dust grains, and therefore any Hα emission from such an accreting protoplanet should be minimally extincted. Indeed, detecting actively accreting protoplanets through Hα emission is the primary motivation behind the GAPlanetS campaign, and this method has been successful twice before (Close et al. 2014; Sallum et al. 2015b).

Certain aggressive KLIP reductions of the 2014 April 12 MagAO data also reveal a point source candidate at the location of HD 100546 "c;" however, a similar point source is also present in the continuum image in all cases, which makes the Hα point source immediately suspect, as we do not expect any significant continuum contribution from a substellar object. Scattered light, on the other hand, should appear the same in Hα and the continuum, and, upon correcting for the Hα excess of the primary star (the source of the light to be scattered), should be fully removed by the SDI process. Indeed, as the SDI-processed images for both data sets reveal, there is no excess in the Hα channel at this location.

In fact, the MagAO images shown in Figure 3 provide an excellent demonstration of the effects of aggressive PSF processing on extended disk structures. There is significantly less rotation in the 2015 May 15 data set than in the 2014 April 12 data set, making the PSF-exclusion criterion necessarily more aggressive (smaller). As a result, structures that appear smooth and extended in the upper panel of the figure appear clumpy and in some cases point-like in the lower panel.

Taken together, these two lines of evidence are consistent with the hypothesis that the source detected at the location of candidate "c" is a scattered-light disk artifact enhanced relative to the disk knot on the opposite side of the major axis by the merger of the eastern inner disk rim with the northeastern spiral arm. Aggressive data processing appears to be the main culprit making this disk feature appear point-like in some reductions.

As a further test of the detectability of an Hα point source in this data, we computed an SDI contrast curve, as shown in Figure 8. To compute this curve, we first convolved the final KLIPed image shown in Figure 3 by a hard-edged circular aperture with a diameter equivalent to the FWHM of the VisAO optical ghost (6 pixels, 004). This ghost serves as an estimate of the unsaturated PSF of the central star and therefore the size of an independently sampled region in the image. The convolved image was divided into annuli with widths equivalent to this measured stellar FWHM. Within each annulus, as many independent apertures as would fit in the annulus without overlapping were placed with a random starting point within the annulus, and the central values in these apertures were recorded. The standard deviation of these central values was taken, multiplied by $1/\sqrt{1+1/n}$ (where n is the number of independent apertures) to account for small sample statistics following Mawet et al. (2014), and multiplied by 5 to generate the 5σ limit for each annulus.

Zoom In Zoom Out Reset image size

This procedure was repeated 500 times (for 500 random realizations of aperture placements) for each annulus, and the values averaged together. To translate this 5σ noise value into contrast, each value was divided by the stellar peak. As HD 100546A was saturated in this data set, the stellar peak was estimated from a measurement of the ghost peak. Using Moffat fits to the stellar and ghost peaks in five unsaturated GAPlanetS data sets, the ghost was shown to have an intensity equivalent to 0.42% ± 0.08% of the stellar peak, and can be scaled by this amount to estimate the stellar peak.

Finally, throughput was computed by injecting fake planets into the raw Hα line images, subtracting the scaled continuum images, and then processing the SDI images with KLIP and the same parameters as the final SDI image. Throughput at a given location is measured as the ratio of the peak brightness of the recovered false planet to the injected planet. The 5σ contrast values were multiplied by this throughput to create the final curve. The curve suggests that we could have detected planets up to ∼1 × 10⁻³ contrast at the location of the HD 100546 "c" candidate and ∼1 × 10⁻⁴ contrast at the location of HD 100546 b.

HD 100546 b is heavily embedded in the disk. Currie et al. (2015) estimate the H-band extinction at the location of the point source candidate to be 3.4 mag, which translates to 22 mag of extinction at R (and therefore Hα) following standard Milky Way extinction laws. This is enough to make any constraints on the accretion luminosity of b meaningless, as we discuss in more detail in the companion to this paper.

The "c" candidate, however, is hypothesized to lie at or near the outer edge of the inner disk rim. If it is heavily embedded in the rim (an unlikely hypothesis given the continuity of disk features at this location), then it suffers from the same problem as b in that dusty material extincts very efficiently at Hα and quickly makes accretion luminosity estimates for embedded protoplanets moot. If the candidate identified by Currie et al. (2015) or hypothesized by Brittain et al. (2014) lies inside the cleared central cavity, however, then the contrast limit at this location can be used to place more meaningful limits on the accretion luminosity and accretion rate of any forming protoplanets, albeit with a number of assumptions as detailed below.

We begin by assuming that the HD 100546 cavity is fully cleared of visible light extincting grains, and indeed the precipitous drop in the Y-pol radial profile approaching the coronagraph supports this assumption somewhat. We take the measured V-band extinction toward HD 100546A (A_V = 0.15; Sartori et al. 2003) and translate it to A_R = 0.11 mag following standard extinction laws (Cox 2000). Following Close et al. (2014), we use this R-band extinction estimate and measured contrast, the zeropoint and width of the Hα filter, and the distance to HD 100546 to translate the measured contrast to an Hα luminosity of $1.57\times {10}^{-4}\,{L}_{\odot }$. If we then assume that empirically derived ${L}_{{\rm{H}}\alpha }$ to L_acc relationships for low-mass T Tauri stars also apply to lower mass objects, then following Rigliaco et al. (2012), this translates to an accretion luminosity of $0.33 \% \,{L}_{\odot }$. Translation of this quantity to an accretion rate requires assumptions about the mass and radius of the accreting object, and we adopt $1.55{R}_{{\rm{J}}}$ and 2M_J in this calculation as reasonably representative of the population of planets we might expect to sculpt the disk rim. Then, following Gullbring et al. (1998), the accretion luminosity translates to an approximate accretion rate of $\dot{M}\approx 1\times {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, corresponding to growth of a Jupiter-mass planet in 100,000 years. The accretion rate onto the primary star is estimated at $\sim {10}^{-7}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Mendigutía et al. 2015), placing our limit at ${\dot{M}}_{\mathrm{planet}}\lt 0.1{\dot{M}}_{\mathrm{star}}$. We note that a number of assumptions have gone into this estimate, including that accretion onto protoplanets happens in a steady flow of material and not stochastically, and thus it is likely only accurate to within 1–2 orders of magnitude.

To examine the effects of our data processing procedures on spiral arms, we produce synthetic images of planet-driven spiral arms in disks using combined hydrodynamics and radiative transfer simulations, and process the simulated images using our GPI pipeline. We adopt the 3 ${M}_{{\rm{J}}}$ planet model of Dong et al. (2016a) with only minor modifications, and briefly summarize salient aspects of the models here. The simulations are described in detail in Dong et al. (2016a; see also Dong et al. 2015, 2016b; Fung & Dong 2015). The simulations are of spiral arms driven by an outer planetary perturber and do not include an inner disk cavity, though we note the qualitative similarity of spiral arms driven by inner and outer planets demonstrated in another work (Zhu et al. 2015). We note that the disk models were adopted without modification, and the location of the planetary perturber does not coincide with the location of the HD 100546 b protoplanet candidate. We leave more precise reproduction of HD 100546's specific disk features, including the inner cavity and prediction of the location of planetary perturbers, to future work.

The three-dimensional density structure of spiral arms in a disk excited by a 3 ${M}_{{\rm{J}}}$ planet was calculated using the code PEnGUIn (Fung 2015). The initial condition of the disk is ${\rm{\Sigma }}\propto 1/r$, and $h/r\propto {r}^{0.25}$, where Σ and h/r are the surface density and aspect ratio in the disk, and h/r at the location of the planet is set to 0.15. The viscosity in the simulation is parametrized using the Shakura & Sunyaev (1973) α prescription with $\alpha =0.01$. The simulation is run for 50 orbits, not long enough for the gap to be fully opened, but sufficiently long for the spiral arms to reach steady state. The resulting 3D disk density structure is subsequently fed into a Monte Carlo radiative transfer code (Whitney et al. 2013) to produce synthetic H-band total intensity images at various inclinations. We convert the gas density as calculated in the hydro simulation to dust density used in the radiative transfer simulation, assuming the dust and the gas are well-mixed, and we adopt the interstellar medium dust model (Kim et al. 1994) for the dust. These dust grains are submicron in size, as assumed in previous scattered light spiral arm modeling works (e.g., MWC 758, Dong et al. 2015; HD 100453, Dong et al. 2016b).

Planet-induced spiral arms are very robust in scattered light imaging independent of the grain properties assumed in the modeling. Qualitative comparisons such as those we are making here are not sensitive to grain models as long as there is small (∼micron-sized) dust present in the disk, as modeling of HD 100546's SED suggests is the case (e.g., Tatulli et al. 2011). Additionally, since small grains dominate the opacity at visible and NIR wavelengths and make up the majority of the dust grains in the surface layers of the disk where scattering originates, the assumption of ISM-like dust properties is reasonable.

To understand the impact of the data processing and to qualitatively assess the reality of features identified around HD 100546, the H-band disk model was convolved with a GPI H-band PSF, injected into a GPI datacube of a disk-free star with comparable brightness and a similar amount of on-sky rotation, and then processed via cADI in precisely the same way as the HD 100546 GPI H-spec data to create a forward model. The underlying surface density model, an H-band total intensity model generated via Monte Carlo radiative transfer modeling as described above, and the forward-modeled image are shown alongside the actual on-sky HD 100546 cADI H-spec image in Figure 9.

Zoom In Zoom Out Reset image size

The forward-modeled image suggests that a two-armed spiral disk perturbed by a single planetary companion and viewed at moderate inclination can result in observed structures that are similar in location, number, brightness, and extent to the features that we observe in HD 100546. This experiment serves as a first-order, albeit striking, demonstration of similarity, and we leave more precise matching and derivation of disk and planet properties from forward models for future work.

The disk models also naturally produce a near-side bright lane feature offset from the rest of the disk and similar in morphology to the feature labeled "A" in Figure 6. Physically, it corresponds to the outer edge of the bottom side (opposite the disk midplane relative to the rest of the disk emission on both near and far sides) of the disk, and the dark region between it and the other disk features corresponds to the dense disk midplane. This bright lane feature in the raw model and forward-modeled images is beyond the edge of the image in Figure 9, but it can be seen clearly in Figure 8 of Dong et al. (2016a). Tunable model parameters like the thickness of the disk midplane and the scale of the spiral arms could conceivably bring the top-side features and the bottom-side bright lane feature closer together in the modeled images, as they appear to be in HD 100546, but we leave this for future work. Alternatively, bright lane "A" may correspond to a different variety of disk features altogether.

Both the forward-modeled and observed images show multiple spiral features, the majority of which lie on the back-scattering far side (NE) of the disk. Self-subtraction is clearly seen breaking single spirals from the raw model image into multiple arcs in the forward model, suggesting that several of the features we identified in Figure 6 may belong to contiguous structures. Thus, the forward model also serves to demonstrate the tendency of aggressive PSF-subtraction techniques to create apparent disk clumps along extended features that are smooth in reality, something that will be very important to account for in future studies of planets embedded in circumstellar disks.