Separating extended disc features from the protoplanet in PDS 70 using VLT/SINFONI

TL;DR

This study uses advanced imaging techniques to analyze the PDS 70 transition disc, successfully redetecting the protoplanet and extended features, and comparing methods for optimal detection and characterization.

Contribution

The paper introduces an optimized PCA-based PSF subtraction method for extended feature recovery and adapts the NEGFC technique for precise astrometry of the protoplanet.

Findings

Redetection of PDS 70 b and disc features

PCA-SADI better recovers extended features

Derived astrometry and spectro-photometry of the protoplanet

Abstract

Transition discs are prime targets to look for protoplanets and study planet-disc interactions. We present VLT/SINFONI observations of PDS~70, a transition disc with a recently claimed embedded protoplanet. We take advantage of the angular and spectral diversity present in our data for an optimal PSF modeling and subtraction using principal component analysis (PCA). We report the redetection of PDS 70 b, both the front and far side of the outer disc edge, and the detection of several extended features in the annular gap. We compare spectral differential imaging applied before (PCA-SADI), and after (PCA-ASDI) angular differential imaging. Our tests suggest that PCA-SADI better recovers extended features, while PCA-ASDI is more sensitive to point sources. We adapted the negative fake companion (NEGFC) technique to infer the astrometry of the companion, and derived $r = 193.5 \pm 4.9 \mathrm{mas}$ and PA = 158.7deg $\pm$ 3.0deg. We used both NEGFC and ANDROMEDA to infer the $K$-band spectro-photometry of the protoplanet, and found results consistent with recent VLT/SPHERE observations, except for their 2018/02 epoch measurement in the $K2$ filter. Finally, we derived an upper limit of $\dot{M_b} < 1.26 \times 10^{-7} \big[ \frac{5 M_{\rm Jup}}{M_b} \big] \big[ \frac{R_b}{R_{\rm Jup}}\big] M_{\rm Jup} $ yr$^{-1}$ for the accretion rate of the companion based on an adaptation of PCA-SADI/PCA-ASDI around the Br$\gamma$ line.

Equations

$$\displaystyle F_{\rm meas}=\big{[}F_{\textrm{Br}\gamma,p}+F_{\textrm{cont},p}+\epsilon_{\textrm{Br}\gamma,\star}\big{]}-\big{[}\frac{F_{\textrm{Br}\gamma,\star}}{F_{\textrm{cont},\star}}F_{\textrm{cont},p}+\epsilon_{\textrm{cont},\star}\big{]}$$

$$\displaystyle F_{\textrm{Br}\gamma,p}=F_{\rm meas}-\bigg{(}1-\frac{F_{\textrm{Br}\gamma,\star}}{F_{\textrm{cont},\star}}\bigg{)}F_{\textrm{cont},p}$$

$$\log (L_{\mathrm{acc}}/L_{\text{sun}})=0.9\times (\log L({\mathrm{Br}}_{\gamma })/L_{\text{sun}}+4)-\mathrm{0.7.}$$

Full text

††pubyear: 2018††pagerange: Separating extended disc features from the protoplanet in PDS 70 using VLT/SINFONI–C

1 Introduction

Transition discs (TDs) are protoplanetary discs showing evidence of inner clearing. These inner clearings were first inferred from the SED, showing a lack of near- to mid-IR excess (e.g. Strom et al., 1989), and recently confirmed by resolved sub-mm observations (e.g. Andrews et al., 2011). Owen (2016) recently highlighted the distinction between TDs with large gaps/cavities ($\gtrsim 20$ au), which show both high accretion rates and bright sub-mm fluxes, and mm-faint TDs which harbour smaller gaps. TDs with large gaps/cavities appear incompatible with a photo-evaporated inner disc, with the cavity instead likely carved by embedded companions. Yet, finding direct evidence of forming planets has proved difficult in these discs despite multiple claims (Kraus & Ireland, 2012; Quanz et al., 2013; Reggiani et al., 2014; Biller et al., 2014; Currie2015; Sallum et al., 2015; Reggiani et al., 2018). The problem is that aggressive filtering on bright and extended disc emission can create point-like artefacts, which could be confused with substellar companions, as seems to be the case in HD 169142 (e.g. Ligi2018), in LkCa 15 (Thalmann et al., 2016, Currie et al. 2019, ApJL in press), and perhaps also in HD 100546 (Follette et al., 2017; Rameau et al., 2017).

While direct observational confirmation of giant planets in transition discs has turned out to be more difficult than expected, there is mounting indirect evidence of their presence (e.g. Casassus, 2016; Huang et al., 2018). Large gaps of different sizes at near-IR and sub-mm wavelengths (e.g. Garufi et al., 2013), asymmetric dust distributions (e.g. Casassus et al., 2013; van der Marel et al., 2013) and spiral arms (e.g. Muto et al., 2012; Christiaens et al., 2014; Benisty et al., 2015) can be produced by a range of different mechanisms, but they are all potential by-products of companion-disc interactions. Price et al. (2018) recently showed that in the case of the transition disc around HD 142527, the characteristics of the disc including the gap size, the banana-shape mm-dust trap and the spiral arms at the gap edge, could all be qualitatively accounted for by the low-mass binary companion (e.g. Biller et al., 2012; Lacour et al., 2016). Could massive companions be at the origin of similar features observed in other transition discs?

High-contrast imaging is one of the most suitable techniques to detect faint signals such as newborn planets or circumstellar disc features in the close vicinity of young stars (Absil & Mawet, 2010; Bowler, 2016). Differential imaging techniques allow to optimize the contrast reached at small angle. Here, we focus on the combination of angular (ADI; Marois et al., 2006) and spectral differential imaging (SDI111It was first referred to as spectral deconvolution (SD) in those works to avoid confusion with dual-band (simultaneous) spectral differential imaging.; Sparks & Ford, 2002; Thatte et al., 2007) applied to medium-spectral resolution integral field spectrograph (IFS) data obtained with VLT/SINFONI on the transition disc of PDS 70.

PDS 70 is a TTauri star surrounded by a disc with a remarkably large dust-depleted inner region. Previous literature has used a distance of $\sim$140 pc, assuming membership of PDS 70 to the Upper Centaurus Lupus association (Riaud et al., 2006). However, the new Gaia DR2 distance places it at only $113\pm 1$ pc (Gaia Collaboration et al., 2018), which alter previous stellar parameter estimates. The near-IR JHK spectrum of the star is compatible with an M1V pre-main sequence star (Long et al., 2018), although UBVRI photometry suggests a K5V star (Gregorio-Hetem & Hetem, 2002). Based on a Markov-Chain Monte-Carlo method, Müller et al. (2018) provided new age and mass estimates of $5.4\pm 1.0$ Myr and of $0.76\pm 0.02M_{\odot}$, resp., for the central star. Considering the new distance throughout this paper, the size of the $\mu$m-size and mm-size dust cavities are $\sim$53 au and $\sim$65 au (Hashimoto et al., 2012, 2015; Long et al., 2018). Near-IR polarimetric images suggest an outer disc inclination of 45-50° and a PA of the semi-major axis of $\sim$159°(Hashimoto et al., 2012). SED modeling hints at the presence of a small optically thick inner disc (Dong et al., 2012). The different gap sizes for $\mu$m-size and mm-size dust are suggestive of the presence of one or more companions inside the large annular gap (e.g. Pinilla et al., 2012; de Juan Ovelar et al., 2013). A point source has been recently detected in the annular gap using VLT/SPHERE data, and interpreted as a forming planet (Keppler et al., 2018; Müller et al., 2018, hereafter K18 and M18, resp.). The point source was tentatively re-detected in the H$\alpha$ line, suggesting on-going accretion and further supporting the protoplanet hypothesis (Wagner et al., 2018). While this work was under review, Keppler et al. (2019) presented new ALMA observations showing new substructures in the disc, some of which matched the NIR observations of M18.

In this paper, we combine spectral and angular differential imaging to reach high contrast in a VLT/SINFONI dataset acquired on PDS 70. We first present our observations and data calibration procedure (Section 2). We then explain our post-processing algorithms and show the resulting final images (Section 3). We report the redetection of a point source at the location of the protoplanet, the redetection of both the front and far side of the outer disc edge, and the detection of several extended features in the annular gap, some of which present similarities to those claimed in M18. In order to better interpret our results, we injected both synthetic spiral arms and fake companions into our dataset and re-processed it with each of our post-processing algorithms. Section 4 is devoted to the characterization of the protoplanet in terms of astrometry, broad-band photometry, spectro-photometry and Br$\gamma$ emission. Finally, we discuss the disc features identified in our images (Section 5). Our main conclusions are summarized in Section 6.

2 Observations and data reduction

We observed PDS 70 with the AO-fed IFS VLT/SINFONI (Eisenhauer et al., 2003; Bonnet et al., 2004) on May 10th 2014, as part of program 093.C-0526 (PI: S. Casassus). A total of 116 datacubes (NEXP) were obtained, consisting of single exposures (NDIT=1) of 60s (DIT) each, hence amounting to 116min integration. We took advantage of the pupil-tracking mode, which resulted in a total parallactic angle variation of 99.8°. Observations were acquired in clear conditions with an average seeing (0$\aas@@fstack{\prime\prime}$7-0$\aas@@fstack{\prime\prime}$9 at $\lambda=$500 nm). The airmass spanned 1.06–1.20. We used the $H$+$K$ grating (spectral resolution $\sim$1500) and the 12.5mas plate scale, resulting in a 0$\aas@@fstack{\prime\prime}$8 x 0$\aas@@fstack{\prime\prime}$8 field of view in each datacube. However, in order to artificially increase our field of view and better sample the large annular gap in the disc, we followed a four-point dithering pattern throughout the observation, placing the star close to a different corner of the detector in consecutive exposures.

Data reduction followed the same procedure as described in Christiaens et al. (2018) for the HD 142527 dataset obtained the same night with the same observing strategy. We refer to that work for details on SINFONI data calibration. The only difference is that we median-combined each set of 4 consecutive cubes of PDS 70 after placing the star at the center of larger frames, such that the final calibrated cube has dimensions of 2000 x 29 x 101 x 101 (wavelength, time, $y$ and $x$ axes, resp.). For the centering, the centroid of the star was found in each frame by fitting a 2D Moffat function implemented in the Vortex Imaging Pipeline222Available at: https://github.com/vortex-exoplanet/VIP. (VIP; Gomez Gonzalez et al., 2017, Gomez Gonzalez et al. 2019, in prep.).

3 Post-processing

This section focuses on the PCA-based algorithm we implemented to take advantage of both the spectral and angular diversity in our SINFONI data to produce high-contrast final images of the environment of PDS 70. We also used the ANDROMEDA algorithm, which leverages the angular diversity in each individual spectral channels. We defer the interested reader to Cantalloube et al. (2015) for a detailed description of ANDROMEDA.

3.1 Algorithm description

Calibrated frames obtained with SINFONI are affected by more speckles in the stellar halo than images obtained with extreme-AO instruments. Nonetheless, the advantage of SINFONI is its large wavelength coverage in the $H$+$K$ mode ($\sim$ 1.45-2.45 $\mu$m), where the stellar diffraction pattern (including speckles) moves radially by an amount proportional to the wavelength in different spectral channels. To determine the optimal post-processing method for our 4D cube, we tested several variations of the PCA algorithms implemented in VIP, which are themselves based on the original PCA-ADI (Soummer et al., 2012; Amara & Quanz, 2012) and sPCA algorithms (Absil et al., 2013).

We first compared (1) PCA in a single step considering both the angular and spectral variation together, (2) PCA-SDI + PCA-ADI in two consecutive steps (hereafter PCA-SADI), and (3) PCA-ADI + PCA-SDI in two consecutive steps (hereafter PCA-ASDI). For both PCA-SADI and PCA-ASDI, the PCA-SDI part was performed in full frame, while the PCA-ADI part was performed in concentric annuli of 2-FWHM width where, for each annulus, the PCA library only includes frames where any putative companion would rotate by at least 1 FWHM. This frame selection enabled us to minimize self-subtraction of faint signals of interest (Absil et al., 2013). PCA-SADI and PCA-ASDI consistently recovered the well-characterized western (forward-scattered) edge of the outer disc of PDS 70, and for some reduction parameters (detailed in the next paragraph) redetected either the faint far side of the disc (PCA-SADI) or the protoplanet candidate claimed in K18 and M18 (PCA-ASDI). By contrast, PCA in a single step did not unambigously detected the bright forward-scattered edge of the disc. We thus discarded PCA in a single step for the rest of the analysis.

Next, we tested different minor variations for the PCA-SADI and PCA-ASDI algorithms. We specifically tested:

using $K$-band spectral channels only, which benefit from a significantly better AO correction than the $H$-band channels, instead of using all $H$+$K$ channels; 2. 2.

normalizing spectral channels based on the stellar flux measured in a 1FWHM-aperture prior to the PCA-SDI part of either PCA-SADI or PCA-ASDI; 3. 3.

using a minimum radial motion for the PCA-SDI part to only include spectral channels that would limit self-subtraction (see also Thatte et al., 2007); 4. 4.

using a maximum threshold in azimuthal motion when building the PCA-ADI library in order to limit the inclusion of temporally decorrelated speckles; 5. 5.

collapsing residual frames (obtained after PCA modeling and subtraction) using either a simple median or a variance-based weighted average (Bottom et al., 2017), for both the PCA-ADI and PCA-SDI parts. For this purpose, we implemented the SDI equivalent of the algorithm presented for ADI residuals in Bottom et al. (2017).

We determined the optimal combination of reduction parameters based on the recovery of known key features of the disc in our post-processed image: the bright forward-scattered edge of the outer disc (West), the dim back-scattered edge of the outer disc (East), and a bright blob $\sim$0$\aas@@fstack{\prime\prime}$19 to the SE of the star (Hashimoto et al., 2012, \al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018). The conclusions from our tests are as follows:

No PCA-SADI and PCA-ASDI reduction using all $H$+$K$ channels managed to recover the faint back-scattered edge of the outer disc, in contrast to PCA-SADI reductions using only $K$-band channels. This suggests that the speckle subtraction during PCA-SDI is optimal when the high spectral correlation of the better AO-corrected $K$-band channels is not diluted in a larger PCA library. The effect of a better speckle subtraction appears to dominate over both the enhanced self-subtraction (due to less radial motion in SDI) and the lower number of accumulated photons (due to the use of only half of all spectral channels). 2. 2.

Normalizing spectral channels in flux led to a better redetection of the far side of the disc in the case of PCA-SADI, but did not improve the redetection of the blob in PCA-ASDI images. Since we will infer the spectro-photometry of the companion in Section 4.1.3 using an adaptation of PCA-ASDI, we only normalize spectral channels in the case of PCA-SADI. 3. 3.

Using a minimum radial motion for the creation of the PCA-SDI library on $H$+$K$ channels implies the inclusion of comparatively more $H$-band spectral channels in the PCA library of $K$-band channels, and vice versa. This leads to a worse speckle subtraction due to worse AO performance in $H$ band than in $K$ band - and hence corresponding decorrelation of speckles. Limiting our analysis to $K$-band channels, we built our PCA-SDI library using only the shortest wavelength channels of that band which kept a high speckle correlation with all other $K$-band channels and did not suffer from significant telluric absorption (spectral channels from 1.93 to 1.95 $\mu$m). The choice of these channels is further justified in Section 4.1.1. 4. 4.

Using a maximum threshold for azimuthal motion to build the PCA-ADI library generally yields a better speckle subtraction (e.g. Wagner et al., 2018). However, with only 29 temporal cubes spanning $\sim$100° rotation, adding this constraint to the minimum azimuthal threshold significantly decreases the number of principal components that can be used for ADI (e.g. $n_{\rm pc}^{\rm ADI}\lesssim 3$ for a 1–3 FWHM azimuthal motion range). We noticed that using 3–10 principal components without a maximum rotation constraint led to a more conspicuous redetection of the companion using PCA-ASDI (Section 3.2). 5. 5.

We did not obtain a significant improvement in the quality of post-processed images using a variance-based weighted average of residuals compared to the median. However, we later noticed that the variance-based weighted average led to a more accurate recovery of spectra of injected fake companions than the median, using the technique presented in Section 4.1.1.

Our final choice of post-processing parameters for PCA-SADI and PCA-ASDI is thus the following: (i) we only consider $K$-band spectral channels, (ii) we normalize spectral channels in flux before PCA-SADI but not for PCA-ASDI, (iii) we set manually the PCA reference library to spectral channels between 1.93 and 1.95 $\mu$m for the SDI part, (iv) we do not set a maximum azimuthal motion to select frames for the ADI part (but do set a minimum threshold equivalent to a 1-FWHM linear motion), and (v) we consider the variance-based weighted average of residuals (for both the ADI and SDI parts) of channels longward of 2.0$\mu$m,

Finally, in addition to PCA-SADI and PCA-ASDI, we also tested PCA-ADI in concentric annuli on the two cubes obtained after collapsing the spectral channels of the $H$ and $K$ bands separately. No significant signal was found in the PCA-ADI $H$-band image, likely due to the poorer AO correction. However, this reduction still enabled us to set an independent upper limit on the $H$ band flux of the companion (Section 4.1.1). In contrast, PCA-ADI on the collapsed $K$ band yielded a tentative redetection of both the forward-scattering edge of the outer disc and the companion.

3.2 Final images

Figures 1a and 1b show the images obtained with our final choice of post-processing parameters for PCA-SADI and PCA-ASDI (Section 3.1), respectively. In order to assess the significance of features in our images, we show standardized trajectory intensity mean (STIM) maps (Pairet et al., 2018a) in the lower two panels of Figure 1, preferring these over classical SNRt maps (Mawet et al., 2014). This choice is detailed in Appendix A. STIM maps are 2D detection maps defined at each location ($i$,$j$) as $\mu(x_{i,j})/\sigma(x_{i,j})$, where $\mu(x_{i,j})$ and $\sigma(x_{i,j})$ are the mean and standard deviation of trajectory $x_{i,j}$ throughout the derotated cube of residual images (i.e. the cube of images obtained after PCA modeling and subtraction, and subsequent alignment with North up and East left), respectively. Using STIM maps overcomes the difficulties of SNRt maps to identify significant signals at short separation from the star and in presence of bright disc features spanning a range of radii (as they bias the noise level estimation at each radius).

The STIM maps corresponding to the final PCA-SADI and PCA-ASDI reductions are shown in figure 1c and d, respectively. The maximum value of the STIM map color scales corresponds to the maximum pixel value obtained in the respective inverse STIM maps (Figure 11). The inverse STIM map is the detection map obtained by derotating the residual images using opposite parallactic angles. Since this procedure does not sum constructively authentic circumstellar signal while preserving the time dependence of residual speckles, it enables us to estimate the maximum pixel value that can be obtained in the detection map purely as a result of combining residual speckle noise (Marois et al., 2008; Wahhaj et al., 2013). In other words, pixels with values above that threshold in our (correctly derotated) STIM maps are unlikely to trace spurious signals. The threshold values are 1.1 and 1.8 in the PCA-SADI and PCA-ASDI STIM maps, respectively. In Figure 1c and d, we identify Features i, ii and iii and b as significant, i.e. containing pixels above that threshold. We identify additional tentative Features iv, v, and vi, (i.e. containing pixels slightly below that threshold) based on the presence of similar features in the long-integration $K$-band image obtained with SPHERE/IRDIS (M18).

Feature b coincides with the protoplanet claimed in \al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018 and Wagner et al. (2018). Feature i corresponds to the bright forward-scattered edge of the outer disc, which has been identified in previous observations of the disc (e.g. Hashimoto et al., 2012; Long et al., 2018, K18). Feature ii has been identified as a possible gap-crossing bridge in M18 (a blend of Features (2) and (3) in their figure B1). Feature iii appears to connect the outer disc and the protoplanet, but has no correspondence in the images of M18. Features iv and v appear to match signals in the vicinity of Features (3) and (5) in figure B1 of M18. Finally, Feature vi is consistent with the dim back-scattered light from the far side of the outer disc, also reported in M18. Feature b is analyzed in details in Section 4, while extended features are further discussed in Section 5.

In Figure 1, the choice of the number of principal components ($n_{\rm pc}$) is such that the detection of the faint back-scattered edge of the outer disc and the redetection of the protoplanet candidate are optimal in the PCA-SADI and PCA-ASDI images, respectively. For comparison, images obtained with other values of $n_{\rm pc}^{\rm ADI}$ are shown in Figure 2. All the features identified in Figure 1c and d are also qualitatively recovered for all values of $n_{\rm pc}^{\rm ADI}$ from 1 to 10, albeit with increasing self-subtraction affecting extended features for increasing $n_{\rm pc}^{\rm ADI}$. In particular, the faint far side of the disc can only be seen for the lowest $n_{\rm pc}^{\rm ADI}$ values for PCA-SADI, but is hardly seen in all PCA-ASDI reductions. PCA-ASDI images with low $n_{\rm pc}^{\rm ADI}$ values also appear to leave stronger residuals at the edge of the 0$\aas@@fstack{\prime\prime}$16 numerical mask. Only PCA-ASDI reductions with sufficiently large values of $n_{\rm pc}^{\rm ADI}$ ($\geq$ 3) enable to redetect the protoplanet candidate claimed in \al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018 and Wagner et al. (2018). Nonetheless, these images also show some extended structures at $\sim$ 0$\aas@@fstack{\prime\prime}$2 to the NNE and at$\sim$ 0$\aas@@fstack{\prime\prime}$3 to the SE of the image, none of which being recovered in the PCA-SADI images.

Although not shown here, the results obtained with $n_{\rm pc}^{\rm SDI}=2$ are very similar to those obtained with $n_{\rm pc}^{\rm SDI}=1$. Increasing $n_{\rm pc}^{\rm SDI}$ to larger values leads to increasing self-subtraction.

Altogether, Figures 1 and 2 suggest that PCA-SADI with minimal number of principal components yields the best subtraction of speckles while minimizing self-subtraction of authentic signals, as it recovers better the faint far side of the disc. Nonetheless, the fact that only PCA-ASDI with large $n_{\rm pc}^{\rm ADI}$ recovers a bright blob at the location of the previously claimed protoplanet candidate is puzzling. Could Feature $b$ trace (at least partially) a co-located extended disc feature, possibly filtered into a point-like source by ADI processing? To help answer this question, we performed two sets of tests; post-processing our data after injection of either extended or point-like features in our pre-processed data cube. The results of these tests are presented in the next two sections.

3.3 Effect of post-processing on spiral arms

3.3.1 Geometric biases

To evaluate how our PCA algorithms filter out extended signals, we injected fake spiral arms in copies of our calibrated 4D cube and re-processed these cubes with PCA-SADI, PCA-ASDI and PCA-ADI (on the collapsed $K$-band spectral channels). Results are shown in Figure 3. We base the shape of the injected synthetic spirals on the best-fit of the trace of Feature iii (Figure 1c) to the linear spiral density wave equation, obtained when forcing the perturber location to be that of protoplanet candidate PDS 70 b (e.g. Rafikov, 2002; Muto et al., 2012). The exact shape of the two injected synthetic spiral arms correspond to that spiral model flipped with respect to the x-axis (spiral #1; Figure 3d), and rotated by 45° clockwise (spiral #2; Figure 3h), respectively. The synthetic spirals are injected at roughly the same contrast level as Feature iii. We considered a flat contrast at all wavelengths (i.e. the spectrum of the injected spirals are a scaled down version of the stellar spectrum). For comparison, the first row of Figure 3 corresponds to the images obtained without fake spiral injection.

The PCA-SADI images with $n_{\rm pc}^{\rm ADI}=1$ and $n_{\rm pc}^{\rm SDI}=1$ recover visually the full trace of the two injected synthetic spirals (Figure 3e and i). By contrast, PCA-ASDI filters out parts of the injected spirals. These are indicated by cyan arrows in Figure 3f and j. In particular, for both spirals, the cyan arrows appear roughly in the middle of the portion of the spiral where the pitch angle is roughly zero (i.e. where the arm is most similar to a circle). We also tested other values of $n_{\rm pc}^{\rm ADI}$ than shown in Figure 3 and noticed that PCA-ASDI filters out spiral #1 at the location of the cyan arrow in Figure 3f for all values of $n_{\rm pc}^{\rm ADI}$ from 1 to 10, while spiral #2 is filtered out at the location of the cyan arrow in 3j for $n_{\rm pc}^{\rm ADI}\geq 3$. By contrast, for all values of $n_{\rm pc}^{\rm ADI}$ from 1 to 10, PCA-SADI filters out neither spiral #1 nor spiral #2. The only effect of increasing $n_{\rm pc}^{\rm ADI}$ in PCA-SADI is to enhance self-subtraction relatively homogeneously along the spiral trace, so that the full shape is preserved. Compared to PCA-SADI and PCA-ASDI, PCA-ADI alone does not appear to optimally subtract speckles, as only parts of the synthetic spirals coincident with a bright speckle or the forward-scattered disc edge produce enhanced blobs in the final image (see Figure  3g).

The results of this test are in agreement with the survival of extended features identified in our final PCA-SADI image (Figure 1a and c) when increasing $n_{\rm pc}^{\rm ADI}$, while not being recovered in the PCA-ASDI images (Figure 2). We also noticed that PCA-SADI with $n_{\rm pc}^{\rm ADI}=1$, $n_{\rm pc}^{\rm SDI}=1$ does not appear to introduce significant geometric biases, as the shape of the retrieved spirals is very similar to that of the injected ones. This suggests the extended features identified in our final PCA-SADI image are real.

3.3.2 Surface brightness and throughput

Next, we quantify the amount of signal that is preserved by our algorithms (i.e. the throughput) along each injected spiral. The latter are injected at a constant surface brightness of 9.84 mag/arcsec2 (77.3 mJy/arcsec2 at $K$ band) throughout their full trace (Figure 3d and h). Figure 4 shows the throughput profile measured in the final PCA-SADI and PCA-ASDI images (blue and red curves, resp.) along the spirals. Since the injected spirals overlap with authentic extended disc features, we subtracted the flux measured along the spiral traces in the images obtained without spiral injection (i.e. Figure 3a and b for PCA-SADI and PCA-ASDI resp.) to that measured in the images obtained after spiral injection (Figure 3e, f, i and j) before computing the throughput.

We notice that both algorithms lead to low throughput. Values $\lesssim$50% are obtained at all radial separations smaller than $\sim 0\aas@@fstack{\prime\prime}3$, with a decreasing trend from large to small radii. This is not surprising: our wavelength range choice (1.95–2.45$\mu$m) leads to insufficient radial motion during the PCA-SDI part of the processing such that partial self-subtraction is unavoidable. Nonetheless, PCA-SADI appears to consistently lead to a higher throughput at all radial separations than PCA-ASDI. Furthermore, the PCA-ASDI profiles show a significant dip at $\sim 0\aas@@fstack{\prime\prime}28$ radius, which is consistent with the dips that were visually identified in Figure 3f and j (cyan arrows).

The spiral injection test also allows us to infer the sensitivity achieved by PCA-SADI towards extended structures in our dataset. The injected spirals are indeed rotated/mirror versions of Feature iii, which is identified in Figure 1c at roughly the limit corresponding to significant signals. We reach thus a sensitivity of $\sim 77$ mJy/arcsec2 in $K$-band total intensity, at $\sim$0$\aas@@fstack{\prime\prime}$25 separation. In comparison, Hashimoto et al. (2012) report a sensitivity of $\sim$0.2 mJy/arcsec2 in $H$-band polarized intensity beyond $\sim 0\aas@@fstack{\prime\prime}25$, and a polarization fraction of $\sim 0.5$%. The sensitivity of PCA-SADI towards extended features in our SINFONI data is thus only slightly lower to what is obtained with polarimetric differential imaging (PDI; e.g. Kuhn et al., 2001; Quanz et al., 2011), although PDI observations at the same wavelength are required for a better comparison.

Finally, we note that at the tip of spiral #1, where the latter intersects the edge of the outer disc, the surface brightness measured in Figure 3e is twice the value measured in Figure 3a. This suggests that Feature i has a similar surface brightness of $\sim 77$ mJy/arcsec2 as the injected spirals. Assuming that Rayleigh scattering is dominant for very small grains, the cross-section is $\propto\lambda^{-4}$ and hence our $K$-band measurement would lead to a prediction of $\sim(1.67/2.18)^{-4}\times 77\approx 224$ mJy/arcsec2 in $H$-band total intensity. This agrees with the $H$-band polarized intensity of $\sim$ $1\pm 0.2$ mJy/arcsec2 and polarization fraction of $\sim 0.5$% measured in Hashimoto et al. (2012).

3.4 Effect of post-processing on point sources

While the results of the spiral injection tests suggest that PCA-SADI reaches a higher sensitivity for faint extended structures, the fact that only the PCA-ASDI images recover a point-like source at the location of claimed protoplanet candidate PDS 70 b raises the question of whether it corresponds to a filtered extended structure, or whether PCA-ASDI is simply better at recovering point-like sources than PCA-SADI. In order to answer this question, we injected two pairs of fake companions in copies of our original 4D cube (Figure 5d and h), and post-processed them with PCA-SADI, PCA-ASDI and PCA-ADI (on the collapsed $K$-band channels). For both tests, companions were (retrospectively) injected based on the contrast as a function of wavelength inferred for the protoplanet candidate (Section 4.1.3). In Test#1, the contrast of the injected companions was rigorously identical to that inferred for PDS 70 b ($\sim 7.6$mag contrast across the $K$ band), while in Test#2 the fake companions were injected at a twice fainter flux in all channels ($\sim 8.4$mag contrast across the $K$ band).

PCA-ASDI recovered all injected fake companions for all $n_{\rm pc}^{\rm ADI}$ values in Test #1, and for all $n_{\rm pc}^{\rm ADI}\geq 3$ in Test #2 (similar to the $n_{\rm pc}^{\rm ADI}$ values required to recover the protoplanet candidate PDS 70 b). By contrast, PCA-SADI with $n_{\rm pc}^{\rm ADI}=1$ and $n_{\rm pc}^{\rm SDI}=1$ is only able to marginally recover both fake companions of Test #1 (Figure 5e). For increasing values of $n_{\rm pc}^{\rm ADI}$, the detections become more significant, although the achieved SNR values are always lower than obtained with PCA-ASDI. For test #2, PCA-SADI only marginally recovers the NE fake companion, but does not recover the other fake companion, for all values of $n_{\rm pc}^{\rm ADI}$. One of the two brighter fake companions (Test #1) is also recovered by PCA-ADI (Figure 5g), while the other fake companions appear at a similar level as the speckle noise.

ADI is known to produce characteristic negative side lobes alongside point-source detections. A hint of these signatures can be seen in Figure 5g. PCA-ASDI also produces negative side lobes (Figure 5f and j), confirming the importance of the ADI part in that algorithm. The low significance of the fake companions in the PCA-SADI images does not enable identification of such features.

3.5 Discussion

We have explored different post-processing techniques and reduction parameters in order to make best use of the angular and spectral diversity present in our SINFONI data. The application of either SDI or ADI had been previously investigated individually on SINFONI data (Thatte et al., 2007; Meshkat et al., 2015, resp.). Our work is the first attempt to use both.

We interpret the poorer results obtained with PCA in a single step compared to our two-step algorithm (PCA-SADI and PCA-ASDI) as speckles being significantly more correlated spectrally than temporally. Using two steps takes better advantage of the different speckle correlation levels, while a single larger PCA library dilutes the high spectral cross-correlation which enables an efficient speckle subtraction in the PCA-SDI step. PCA-SADI and PCA-ASDI appear thus to be the most appropriate methods for SINFONI datasets. This conclusion might not apply for extreme-AO fed IFS such as those of SPHERE/IFS, GPI or SCExAO/CHARIS because 1) the better AO system provides more consistent high Strehl ratio images as a function of time, and 2) the IFS configuration is different (lenslet-based vs image-slicer based for SINFONI) so that the spectral correlation of speckles might be different (e.g. Claudi et al., 2008; Wolff et al., 2014; Brandt et al., 2017). Applying our two-step PCA algorithms led to the recovery of features identified in recent observations of the system using extreme-AO instrument SPHERE (Figure 1; \al@Keppler2018, Muller2018; \al@Keppler2018, Muller2018, ).

Using PCA-SDI in the first step is expected to better subtract speckles than PCA-ADI because of the stronger spectral correlation of speckles (acquired simultaneously in all different spectral channels) than the temporal correlation throughout the observing sequence. The latter is subject to both slowly varying quasi-static speckles and variations in observing conditions (see also Rameau et al., 2015). This interpretation is consistent with the fact that only 1 principal component for both the SDI and ADI parts ($n_{\rm pc}^{\rm SDI}=n_{\rm pc}^{\rm ADI}=1$) is enough for PCA-SADI to identify all extended structures in the disc, while increasing values of either $n_{\rm pc}^{\rm SDI}$ or $n_{\rm pc}^{\rm ADI}$ progressively self-subtract them while not revealing any new feature. By construction, PCA-SDI preserves azimuthally extended structures, while being potentially more harmful for radially extended structures. However, our final reductions only use the $K$-band spectral channels (instead of the whole $H+K$ range), hence the correspondingly lower amount of radial motion of authentic circumstellar signal (disc or planet) during the SDI part might still induce some (radial) self-subtraction.

Performing PCA-ADI first whitens the noise, hence PCA-SDI in the second step is expected to be less efficient. Whitening is not perfect though, in particular for small values of $n_{\rm pc}^{\rm ADI}$, as residual speckles still show significant correlation between different spectral channels even after PCA-ADI (e.g. Figure 7). As the temporal correlation is not as high as the spectral correlation, a larger $n_{\rm pc}^{\rm ADI}$ is required for PCA-ASDI than for PCA-SADI to build an optimal model of the speckle pattern.

Post-processing of disc images with ADI is known to introduce geometric biases, filtering out parts of the disc (e.g. Milli et al., 2012; Rameau et al., 2017). PCA-ADI with a large $n_{\rm pc}^{\rm ADI}$ makes this effect worse because it includes more signal from the disc in the principal components (e.g. Pairet et al., 2018b). Similar to the results obtained with ADI, our synthetic spiral injection tests show that PCA-ASDI can also filter out significant parts of the spirals (Figure 3f and j). This is consistent with the fact that ADI has a more dominant role in PCA-ASDI than in PCA-SADI. Nonetheless, while ADI can filter out azimuthally extended features, it can also enhance — for the same reason — the detection of azimuthal asymmetries or point sources. This is likely the main reason why only our PCA-ASDI images detect conspicuously PDS 70 b. Our fake companion injections indeed show that the PCA-ASDI algorithm is more sensitive to faint companions than PCA-SADI (Figure 5) — at the expense of subtracting azimuthally extended disc structures.

Rameau et al. (2015) investigated the sensitivity of an algorithm similar to SADI to detect faint companions and found that it strongly depended on the spectral features of the companion. Therefore in our fake companion injection tests, we considered only artificial companions with the same contrast as a function of wavelength as the protoplanet PDS 70 b (as inferred in Section 4.1.3). It is beyond the scope of this paper to test the recovery of companions with other spectra. However, our conclusions are expected to hold, at least qualitatively, for companions emitting significantly more flux at the long wavelength end of the $K$ band than at the short end.

Our final PCA-ASDI and PCA-SADI images, and the tests carried out in Sections 3.3 and 3.4 alone cannot rule out the possibility that Feature $b$ is the result of filtering of an extended disc feature. To confirm the point-source nature of Feature b, we also applied the ANDROMEDA algorithm (Cantalloube et al., 2015, Cantalloube et al. 2019, in prep.). ANDROMEDA is a maximum likelihood matched-filter algorithm that searches for the expected signature of point-like sources. By construction, it is sensitive to authentic point sources, and not to extended disc features filtered by ADI. Figure 6 shows the detection map obtained after combination of the ANDROMEDA detection maps obtained for all individual spectral channels at wavelength $>2\mu$m, based on a weighted average proportional to the square of the SNR at each pixel (Thiébaut et al., 2016). Feature b is the only signal detected above $3\sigma$ in our image, providing an independent confirmation of the point-source nature of Feature b.

4 Characterization of the protoplanet

In addition to ANDROMEDA’s detection map, several other lines of evidence presented in recent studies argue in favor of Feature b being a protoplanet (see more details in Section 4.4). Therefore, we attempted to extract both the exact astrometry and the spectrum of Feature $b$. We followed two different methods: the negative fake companion technique (NEGFC) combined to an adaptation of our PCA-ASDI algorithm (Section 4.1), and ANDROMEDA (Section 4.2).

4.1 Spectro-astrometry inferred using NEGFC

4.1.1 Negative fake companion technique

NEGFC is a specific form of forward modeling used to characterize faint point sources, and first proposed in the context of ADI datasets (e.g. Marois et al., 2010; Lagrange et al., 2010; Soummer et al., 2012). Our adaptation of the NEGFC algorithm is summarized in Figure 7.

We first attempted to extract the exact position and contrast as a function of wavelength of Feature b using NEGFC coupled with PCA-ADI applied to each spectral channel, as performed for HD 142527 B in Christiaens et al. (2018). Taking advantage of the whitening of the noise by ADI, we binned (40 by 40) adjacent spectral channels after PCA-ADI processing (Figure 7a to d), and over the whole $K$-band wavelength range (Figure 7e). However, we noticed that (1) the pixel values of surrounding bright speckles corresponded to about half of the value of pixels at the location of the companion, and that (2) the contrast estimated by NEGFC (PCA-ADI) corresponded to a roughly twice brighter flux than estimated in the $K1$ and $K2$ filters using SPHERE (M18). This suggested that a significant contribution from underlying speckles was indeed biasing the flux estimates.

In order to better subtract the residual speckle in the image, PCA-ADI was followed by PCA-SDI, similar to the PCA-ASDI algorithm described in Section 3.1. PCA-ASDI images (right column of Figure 2 or Figure 7j) are obtained after PCA-SDI model PSFs are subtracted to each ADI-processed spectral images, and all residual spectral frames of the $K$ band (2.0–2.45$\mu$m) are combined using a variance-weighted average. In comparison, images labelled PCA-ASDIλ (Figure 7f to i) are obtained by combining residual adjacent spectral channels 40 by 40. We notice in Figure 7f to j that the residual speckle noise has been attenuated compared to images obtained with PCA-ADI alone, so that the signal from the companion appears now to stand out more conspicuously.

Since the spectrum inferred by NEGFC (PCA-ASDI) can be biased if the companion is also present in the PCA library used for the PCA-SDI part (Maire et al., 2014; Rameau et al., 2015; Galicher et al., 2018), we carefully selected the spectral channels to be included in the PCA-SDI library as those with: (1) no signal departing from speckle noise at the location of the companion, and (2) showing the most correlated residual speckle noise with spectral channels where the emission of the companion is detected. The first condition is met for all spectral channels shortward of $\sim 2.0\mu$m, while the second condition is only met for $K$-band spectral channels, since the AO correction was significantly poorer in the $H$-band channels. In order to minimize the risk of self-subtraction, we selected the shortest wavelength spectral channels of the $K$ band that were not significantly affected by telluric line absorption, and hence limited the PCA library to include all spectral channels between 1.93 and 1.95 $\mu$m. As can be seen in Figure 7k, the residual speckle noise after PCA-ADI in the 1.93–1.95$\mu$m spectral channels correlates with the residual speckle pattern in spectral channels longward of 2.0 $\mu$m (Figure 7a–d), but no significant signal from the companion is detected. This is consistent with the drop in flux shortward of $\sim$2.1$\mu$m in both the spectrum of the companion and corresponding best-fit models presented in M18.

Our strategy is to first use PCA-ASDI to infer the most accurate astrometry of the companion, taking advantage of the enhanced SNR, and then to use PCA-ASDIλ to extract the contrast as a function of wavelength of the companion, at the exact location inferred with PCA-ASDI. For this purpose, we adapted the version of NEGFC that is implemented in VIP for PCA-ADI (Wertz et al., 2017; Gomez Gonzalez et al., 2017) to work with our PCA-ASDIλ and PCA-ASDI implementations.

4.1.2 Broad-band astrometry

In order to obtain the best estimate of the radial separation ($r$) and PA of the companion, we used NEGFC in combination with PCA-ASDI, as it combines all residual spectral frames. The NEGFC optimization was performed using a Nelder-Mead simplex-based algorithm (see more details in Wertz et al., 2017).

Instead of considering all spectral channels of the $K$ band for PCA-ASDI, we only consider channels in the [2.13—2.29]$\mu$m wavelength range, as they show less residual speckle noise in the vicinity of the protoplanet. We noticed indeed that at shorter wavelength than $\sim 2.13\mu$m, the signal of the companion is too faint, while at longer wavelength than $\sim 2.29\mu$m, a bright extended feature radially inward of the companion appears to slightly bias the centroid estimate, shifting the estimate to about 1 pixel shorter radius. This bright extended feature is roughly shaped as an inner spiral arm stemming from the location of the companion, and is best seen in Figure 7h and i.

Using NEGFC (PCA-ASDI), we find the radial separation and PA of the companion to be: $193.5\pm 4.9$ mas and $158.7\degr\pm 3.0$°. Given the PA of the outer disc (PA$\sim$159°), both the projected and deprojected physical separations are 20.9$\pm$0.6 au. The quoted astrometric uncertainties reflect (1) the error associated with stellar centering; (2) instrumental systematic uncertainties (plate scale, true North and pupil offset); and (3) residual speckle noise. Comparison between centering with 2D Gaussian and 2D Moffat functions are consistent within $\sim$0.05 pixel, and we hence expect the error associated to stellar centering to not be significantly larger than that. We conservatively consider the same systematic uncertainties as estimated in Meshkat et al. (2015) for SINFONI data obtained at a similar epoch, but for a companion at larger separation from its star than PDS 70 b: 0.4 mas and 0.5°. The error associated to residual speckle noise dominate the budget: 4.8 mas and 2.9°, for $r$ and PA respectively. The procedure used to estimate them is similar to that presented in Wertz et al. (2017) and is detailed in appendix B.

4.1.3 Spectro-photometry

We applied NEGFC (PCA-ASDIλ) to determine the contrast of the companion as a function of wavelength. The position of the companion was fixed to that inferred in Section 4.1.2 using NEGFC (PCA-ASDI), leaving only a single free parameter to explore. A grid search was performed to infer the optimal contrast of the companion in each binned spectral channel. Given that the star was not saturated in our datacube, the centering of the frames based on fitting the stellar centroid with a 2D Moffat function is expected to be very accurate. Therefore, the position of the companion is not expected to vary significantly in our post-processed spectral channels.

We did not use the Nelder-Mead simplex-based algorithm on the 3 parameters together in each (binned) spectral channel for the same reasons we limited NEGFC (PCA-ASDI) to spectral channels between 2.13–2.29 $\mu$m for the astrometric estimate in Section 4.1.2. We noticed that for some of the spectral channels outside of that wavelength range, the Nelder-Mead algorithm converged at erroneous locations, hence biasing the inferred spectro-photometry of the companion. The grid search is also faster, which facilitated the estimation of both the residual speckle uncertainties and the optimal $n_{\rm pc}^{\rm ADI}$ used for contrast estimation in each binned spectral channel (appendix B). The final inferred contrast($\lambda$) is shown with blue error bars in Figure 8.

4.2 Spectro-photometry inferred using ANDROMEDA

Based on individual ANDROMEDA detection maps, we selected the spectral channels where the significance of the protoplanetary signal was larger than 3$\sigma$, and built a spectrum based on the flux estimated by ANDROMEDA. These spectral channels were all located in the $K$ band. As we noticed the presence of some potential outliers, we median-combined the fluxes inferred by ANDROMEDA in adjacent channels 20 by 20. The resulting spectrum is shown with yellow error bars in Figure 8. It is consistent with both the spectrum inferred using NEGFC (Section 4.1.3), and the SPHERE/IRDIS measurements acquired in the $K1$ filter ($\sim 2.10\mu$m) on 2016/05/14 and 2018/02/24, but only the $K2$ measurement ($\sim 2.25\mu$m) acquired at the first epoch (\al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018). Most fluxes estimated by either NEGFC or ANDROMEDA between 2.19 and 2.31 $\mu$m (all but two points) are slightly larger than the 2018/02 SPHERE measurement. A detailed spectral analysis of the protoplanet is beyond the scope of this paper and is deferred to a forthcoming paper (Christiaens et al. 2019, submitted to ApJL).

We note that ANDROMEDA performs better than NEGFC+PCA-ASDI for the extraction of the spectrum of the companion. It leads to smaller uncertainties and preserves a finer spectral resolution. This is because ANDROMEDA is not sensitive to extended disc features, while PCA-ASDI does not entirely filter out extended signals from the disc, which complicates the estimation of the flux stemming from the point source only.

4.3 Constraints on Br$\gamma$ emission

Wagner et al. (2018) recently reported the tentative detection of H$\alpha$ emission originating from the location of the protoplanet candidate. In this section, we investigate the presence of Br$\gamma$ line emission (2.16612 $\mu$m) at the location of the companion. We applied two slightly different algorithms for detection of Br$\gamma$ emission, referred to as SADI+Br$\gamma$ and ASDI+Br$\gamma$, and that we adapted from our PCA-SADI and PCA-ASDI algorithms respectively. The two techniques are described in details in Appendix C.

Neither the final ASDI+Br$\gamma$ image nor the final SADI+Br$\gamma$ image (Figure 9c and f, resp.) reveal any significant Br$\gamma$ signal, neither at the location of the companion nor elsewhere in the image. We note that the throughput-corrected residual noise is slightly lower in the final $I_{\textrm{SADI+Br}\gamma}$ image than in the $I_{\textrm{ASDI+Br}\gamma}$ image. We hence focus on the former one in the rest of this section to estimate an upper limit on the Br$\gamma$-line flux of protoplanet candidate PDS 70 b. The hypothetical flux that would be measured in $I_{\textrm{SADI+Br}\gamma}$ at the location of the protoplanet candidate can be expressed as follows:

[TABLE]

where $F_{\textrm{Br}\gamma,p}$ and $F_{\textrm{cont},p}$ are the Br$\gamma$ line and continuum fluxes of the planet, $F_{\textrm{Br}\gamma,\star}$ and $F_{\textrm{cont},\star}$ are the Br$\gamma$ line and continuum fluxes of the star, and $\epsilon_{\textrm{Br}\gamma,\star}$ and $\epsilon_{\textrm{cont},\star}$ are the contribution of residual speckle signal at the location of the planet in Br$\gamma$ line and continuum channels, respectively. Fraction $F_{\textrm{Br}\gamma,\star}$/$F_{\textrm{cont},\star}$ stems from the rescaling of continuum channels performed at the beginning of the SADI+Br$\gamma$ procedure.

The residual speckle noise in $I_{\textrm{Br}\gamma}$ and $I_{\textrm{cont}}$ appears well correlated (Figure 9d and e), which is also confirmed by the absence of residual signal at the spatial frequency expected for speckles (size $\approx$ FWHM) in $I_{\textrm{SADI+Br}\gamma}$ (Figure 9f). Therefore we can consider $\epsilon_{\textrm{Br}\gamma,\star}$ and $\epsilon_{\textrm{cont},\star}$ in Eq. 1) to efficiently cancel out. Rewriting Eq. 1 in terms of Br$\gamma$ line flux of the protoplanet candidate yields:

[TABLE]

Given the absence of significant signal at the location of the companion in Figure 9f, we replace $F_{\rm meas}$ by an upper limit corresponding to 5 times the noise level at the radial separation of the companion. This noise level is estimated from the standard deviation of the fluxes measured in independent apertures at that radial separation, corrected by a factor accounting for small sample statistics (Mawet et al., 2014). Given the efficient speckle removal and whitening of the noise, we consider that 5 times this noise level is a conservative estimate of the 5$\sigma$ contrast reached at the separation of the companion. It corresponds to $\sim 6.3\times 10^{-4}$ ($\sim$ 8.0 mag) contrast. An estimation for $F_{\textrm{cont},p}$ can be obtained from the photometric point directly longward to the Br$\gamma$ line in Figure 8: $\sim$7.6 mag contrast. Since $F_{\textrm{Br}\gamma,\star}/F_{\textrm{cont},\star}\approx 0.99$, the second term of Eq. 2 is thus almost two orders of magnitude smaller than our upper limit on $F_{\rm meas}$, and can be neglected.

Our final upper limit on $F_{\textrm{Br}\gamma,p}$ is $\sim 8.3\times 10^{-20}$ W m*-2*. It is obtained conservatively considering the 5$\sigma$ contrast achieved in our final $I_{\textrm{SADI+Br}\gamma}$ image, the calibrated spectrum of the star reported in Long et al. (2018) and a 7.5Å-FWHM linewidth. Considering a distance of 113 pc, this translates into a Br$\gamma$ luminosity $\log(L(\mathrm{Br}\gamma/L_{\odot})<-7.48$. Calvet et al. (2004) inferred the following empirical relationship to convert $L(\mathrm{Br}\gamma)$ into an accretion luminosity:

[TABLE]

Although this expression was inferred for T-Tauri stars, we apply it to PDS 70 b due to the lack of similar relationship for lower mass objects. We find $\log(L_{\mathrm{acc}}/L_{\sun})<-3.83$. We convert this accretion luminosity to a mass accretion rate using $\dot{M_{b}}=1.25L_{\mathrm{acc}}R_{b}/GM_{b}$ (Gullbring et al., 1998), and find: \dot{M_{b}}<1.26\times 10^{-7}\big{[}\frac{5M_{\rm Jup}}{M_{b}}\big{]}\big{[}\frac{R_{b}}{R_{\rm Jup}}\big{]}M_{\rm Jup} yr*-1*.

4.4 Discussion

Previous claims of protoplanet PDS 70 b detections all made use of ADI (\al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018; Wagner et al., 2018). Considering our PCA-SADI image and our fake spiral injection tests alone, it is not possible to discard the possibility that Feature b is tracing an extended disc feature that is filtered by ADI. Given the past investigations on the effect of ADI on discs (e.g. Esposito et al., 2014; Follette et al., 2017), including the blobs created along the semi-major axis of thin inclined rings (Milli et al., 2012)), our results would actually suggest extreme caution regarding the blob located along the semi-major axis of the transition disc of PDS 70. The location of Feature b at the tip of a spiral-like feature (Feature iii) seen in the PCA-SADI image is also reminiscent of HD 100546, for which Rameau et al. (2017) showed that a spiral seen in RDI images could produce a point source similar to protoplanet candidate HD 100546 b after ADI filtering (Quanz et al., 2013).

However, in contrast to other protoplanet candidates, several lines of evidence argue in favor of the protoplanet hypothesis in the case of PDS 70 b, apart from the point source detection:

the proper motion of the protoplanet is consistent with an object on a Keplerian orbit (M18); 2. 2.

the inferred spectrum of the protoplanet is compatible with synthetic spectral models of young substellar objects (M18); 3. 3.

H$\alpha$ emission is tentatively detected at the protoplanet location (Wagner et al., 2018); and

These arguments motivated us to extract the exact astrometry and spectro-photometry of Feature b (Section 4.1.1). We found the position of the companion ($r=193.5\pm 4.9$mas and PA $=158.7\degr\pm 3.0$°) to be consistent with detections at other epochs (\al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018). The contrast as a function of wavelength is qualitatively consistent with the two photometric points inferred with SPHERE/IRDIS in the $K1$ and $K2$ filters. Our analysis thus further strengthens points (i) and (ii) listed above. A detailed spectral analysis of the companion using both our SINFONI spectrum and the rest of the SED presented in M18 is beyond the scope of this paper, and will be presented in a forthcoming paper.

An alternative method to spectrally characterize the companion is to cross-correlate spectral templates (e.g. Konopacky et al., 2013; Hoeijmakers et al., 2018). The application of this promising technique is beyond the scope of this work, but could provide further spectral information on the protoplanet.

Based on the non-detection of Br$\gamma$ emission at the location of the companion, we constrained its accretion rate to be \dot{M_{b}}<1.26\times 10^{-7}\big{[}\frac{5M_{\rm Jup}}{M_{b}}\big{]}\big{[}\frac{R_{b}}{R_{\rm Jup}}\big{]}M_{\rm Jup} yr*-1*. A major uncertainty in this value is whether the empirical relationship to convert Br$\gamma$ luminosity to accretion luminosity inferred for T-Tauri stars is still valid in the planetary mass regime (Calvet et al., 2004). Regardless, our upper limit is compatible with the mass accretion rate inferred from the tentative detection of H$\alpha$ emission and using a similar T-Tauri-based conversion relationship: $10^{-8\pm 1}M_{\rm Jup}$ yr*-1* for a 5–9 $M_{\rm Jup}$ planet (Wagner et al., 2018).

The uncertainty on the accretion rate measured in Wagner et al. (2018) stems partially from the unknown value of the extinction in the line of sight of the companion. For $A_{V}\sim 3$mag, the measured H$\alpha$ flux of the protoplanet candidate would imply $\dot{M_{b}}\sim 10^{-7.0}M_{\rm Jup}$ yr*-1* (for a 5–9 $M_{\rm Jup}$). Therefore, our upper limit on the accretion rate based on the Br$\gamma$ line – significantly less affected by extinction – suggests that the line of sight towards the protoplanet candidate has a visual extinction $A_{V}<3$mag.

5 Characterization of disc features

Based on the STIM maps of our PCA-SADI and PCA-ASDI images (bottom row of Figure 1), we identified similar features as in the image obtained from a long integration of PDS 70 using extreme-AO instrument VLT/SPHERE (M18): the bright forward-scattered edge of the disc (Feature i), several azimuthally extended features (Features ii, iv and v), the faint back-scattered light from the far side of the disc (Feature vi), and a bright blob (Feature $b$), interpreted as a protoplanet. The root of features ii and iv is consistent with the location of the spur recently observed in dust continuum and 12CO emission (Keppler et al., 2019). The spur could be tracing the tip of a gap-crossing stream. Such stream would physically require the presence of an additional companion, which could account for spiral-shaped features ii/iv seen in our images.

Feature iii is the only feature with no clear counterpart in the M18 image. It is best seen in our final PCA-SADI image, although a hint of its presence can also be seen in the PCA-ASDI images. It appears to stem from the location of the protoplanet and could either trace an outer spiral or a gap-crossing stream connected to the outer disc. We estimate the deprojected pitch angle of Feature iii to be $27\pm 4$° based on the best-fit logarithmic spiral, if it lies in the disc plane. If authentic, why was it not seen in the long-integration SPHERE/IRDIS observation (M18)? The M18 image was obtained using ADI, which is known to create negative azimuthal side lobes around bright features. In our synthetic spiral injection tests (Figure 3), we noticed that a spiral with the same shape as Feature iii was cancelled by PCA-ASDI in the vicinity of the bright disc edge, where the spiral’s pitch angle was also the lowest (close to an arc circle). This could be due to either self-subtraction and/or hiding by the negative lobe of the bright disc edge. As pointed out in Section 3.5, ADI has a dominant role in our PCA-ASDI algorithm. It is thus possible that the feature was cancelled out in the M18 image for the same reason.

Apart from the forward-scattered edge of the outer disc, neither the extended features identified in out SINFONI mages, nor those identified in the median-ADI image of M18, have a clear counterpart in polarized light images (K18). Polarimetric differential imaging (PDI; e.g. Kuhn et al., 2001; Quanz et al., 2011) is often considered the best method to image extended disc features. The absence of extended features identified with ADI or SADI in PDI images of PDS 70 suggests caution. Alternatively, non-optimal observing conditions (mediocre seeing), a shorter integration time and the fact that polarized intensity contains only a fraction of the total intensity, might account for the lack of a conspicuous counterpart in polarized light to the features identified in our images (K18).

Based on previous hydrodynamical and radiative transfer simulations, an embedded giant planet of several Jupiter masses is expected to launch significant spiral density waves in the disc, with a potentially observable arm outward from the planet location (e.g. Zhu et al., 2015; Dong et al., 2015). In particular, simulations of Zhu et al. (2015) show that for discs with large inclinations the spiral is expected to be significantly brighter in full intensity than polarized intensity (see their Figure 17), which would be consistent with the absence of counterpart to Feature iii in polarized light. Alternatively, Feature iii might be tracing thermal emission from shocks in the spiral rather than scattered light (e.g. Lyra et al., 2016), which would also be consistent with the apparently large observed pitch angle ($\sim 27\pm 4$° after disc deprojection). Yet another possibility is that Feature iii is a gap-crossing stream, which can be produced under certain companion mass ratio and orbital eccentricity conditions (e.g. Price et al., 2018). Either the spiral arm or the gap-crossing stream scenario would be compatible with an accreting protoplanet at the location of Feature b (Wagner et al., 2018). Finally, we note that Feature iii roughly follows the inner edge of CO J=3-2 bulk emission in the South part of the disc, as seen with ALMA (Long et al., 2018). However, new ALMA observations at higher angular resolution are required to confirm the possible connections between features seen in NIR and sub-mm wavelength images.

Altogether we cannot exclude the possibility that the extended features seen in our PCA-SADI image are authentic. New observations of PDS 70 making use of reference star differential imaging (RDI; e.g. Lafrenière et al., 2009; Ruane et al., 2017; Rameau et al., 2017) might confirm and better constrain these extended signals. Whether a spiral arm or a gap-crossing stream is expected in the case of PDS 70 needs to be tested with dedicated hydrodynamical simulations, followed with radiative transfer predictions and post-processing.

6 Summary

We presented the first application of combined spectral and angular differential imaging to VLT/SINFONI data. We detailed the techniques we used to leverage both the spectral and angular diversity for an optimal stellar halo modeling and speckle subtraction using principal component analysis (PCA) in the case of a SINFONI dataset on PDS 70. We then compared the results obtained when applying spectral differential imaging before angular differential imaging (PCA-SADI), or vice-versa (PCA-ASDI). STIM maps (Pairet et al., 2018a) were used to identify significant features in our final images. Our PCA-ASDI image reveals a point source at roughly the same location as the recently claimed protoplanet, while our PCA-SADI image shows both azimuthally extended features at a similar radius as the companion, and the faint back-scattered light from the far side of the disc. Most features have counterparts in the SPHERE image, but a possible spiral arm or gap-crossing stream which appears connected to the location of the companion.

In order to better interpret the observed features and disentangle possible geometric biases inherent to the PCA-SADI and PCA-ASDI algorithms, we carried out a series of tests consisting in the post-processing of copies of our original datacube in which we injected either synthetic spirals (Section 3.3) or fake companions (Section 3.4). We concluded that PCA-SADI recovers better azimuthally extended features, while PCA-ASDI is more sensitive to the detection of faint companions (or at least of companions with a similar spectrum as PDS 70 b), hence accounting for the differences in the PCA-SADI and PCA-ASDI images. Therefore, we recommend the use of both PCA-SADI and PCA-ASDI for the post-processing of similar IFS datasets in the future. We also encourage the systematic use of dedicated tests (such as injection of synthetic extended features) when ADI is used.

Based on our PCA-SADI/PCA-ASDI images and spiral injection tests alone, we cannot rule out the possibility that part of the signal from the location of the companion stems from a filtered extended structure. However, based on the independent ANDROMEDA detection at the location of PDS 70 b, we confirm that at least part of the signal comes from a point source. Furthermore, both the spectrum and multi-epoch astrometry presented in M18 argue in favor of an authentic protoplanet. We hence adapted both the negative fake companion technique and our PCA-ASDI algorithm to infer reliable estimates of the astrometry, and found the protoplanet to be at $r=193.5\pm 4.9$mas ($21.9\pm 0.6$ au) and PA $=158.7\degr\pm 3.0$°in our images, consistent with the astrometry reported in K18 and M18 at other epochs. We then extracted the spectro-photometry of the companion using two independent techniques: our negative fake companion algorithm and ANDROMEDA. Both methods led to estimates of the contrast of the companion (with respect to the star) as a function of wavelength that are consistent with each other. However, the ANDROMEDA spectrum boasts both a higher spectral resolution and smaller uncertainties than the spectrum inferred by NEGFC (PCA-ASDIλ), which we attribute to the ability of ANDROMEDA to ignore any bias from signals that are not point-like in nature. The NEGFC and ANDROMEDA spectra both agree with the SPHERE $K1$ photometric estimates obtained on 2016/05 and 2018/02, and the first epoch $K2$ measurement (\al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018). However, the second epoch $K2$ measurement is slightly fainter than our estimates. A detailed spectral analysis will be presented in a forthcoming paper (Christiaens et al. 2019, subm. to ApJL).

Finally, we also adapted our PCA-SADI and PCA-ASDI algorithm to set constraints on the Br$\gamma$ line emission. Assuming a similar Br$\gamma$ luminosity to accretion luminosity relationship as T-Tauri stars, our non-detection of Br$\gamma$ emission sets a limit on the mass accretion rate of the companion: \dot{M_{b}}<1.26\times 10^{-7}\big{[}\frac{5M_{\rm Jup}}{M_{b}}\big{]}\big{[}\frac{R_{b}}{R_{\rm Jup}}\big{]}M_{\rm Jup} yr*-1*. Considering both our non-detection of Br$\gamma$ and the tentative detection of H$\alpha$ emission, we infer that the visual extinction towards the companion should be smaller than $\sim$3.0 mag.

From a technical perspective, we used the full potential of an image-slicer based IFS not originally conceived to reach extremely high contrasts, and obtained final images of similar quality as more dedicated extreme AO-fed instruments such as VLT/SPHERE. Our results suggest that future IFS with a similar design but fed with extreme adaptive optics systems or in space (e.g. JWST/NIRSPEC) may reach even higher contrasts.

Acknowledgements

VC is thankful to Sebastian Perez, Miriam Keppler and Roy van Boekel for constructive discussions. VC also acknowledges Zachary Long and Mike Sitko for sharing their SpeX spectrum of PDS 70, and Alexander Heger for the use of Oxygen. VC and SC acknowledge support from the Millennium Science Initiative (Chilean Ministry of Economy) through grant RC130007. VC and OA acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (ERC Grant Agreement No. 337569) and from the French Community of Belgium through an ARC grant for Concerted Research Action. We also acknowledge funding from the Australian Research Council via DP180104235. This work was supported by the Multi-modal Australian ScienceS Imaging and Visualisation Environment (MASSIVE) (www.massive.org.au).

Appendix A SNR maps versus STIM maps

Despite a widespread use of advanced differential imaging techniques in the HCI community, assessing the significance of signals found in post-processed images remains challenging and depends on assumptions regarding the statistics of residual speckle intensities (e.g. Marois et al., 2008; Mawet et al., 2014; Jensen-Clem et al., 2018). Mawet et al. (2014) suggested to compute the signal-to-noise ratio (SNR) at each pixel as the ratio between the flux measured in a FWHM-wide aperture at the pixel location and the standard deviation of the fluxes measured in independent FWHM-wide apertures at the same radial separation, corrected by a radius-dependent factor reflecting small sample statistics (Student $t$ statistics). SNRt maps computed in such a way are now routinely used. However, they assume the presence of at most one point source at each radial separation from the star. In the case of PDS 70, the bright, well-characterized, edge of the outer disc extends over a continuous range of projected radii from the NW to the S of the star (Hashimoto et al., 2012, \al@Keppler2018,Muller2018; \al@Keppler2018,Muller2018). Furthermore, some additional extended features might also stem from within the large annular gap (M18). These bright extended features affect the flux measurements in the independent apertures, and hence artificially increase the noise level estimated at each radius. As a consequence, according to the classical SNRt maps corresponding to Figure 1a and b (Figure 10a and b, resp.), no significant signal (SNR$>$ 5) is present in our images, despite the presence of visually conspicuous signals coinciding with features previously identified in different datasets.

Therefore, we adapted the classical SNRt definition in order to estimate a more appropriate SNR map in the presence of extended disc features in the image. Instead of estimating the noise at each radius using independent apertures in the final image, we consider the same apertures but in the image obtained using the opposite values of parallactic angles for derotation of the residual frames. The images produced in such a way preserve the time dependence of residual speckle, while not constructively adding authentic circumstellar signal (point-like or extended), and are hence good proxies of the residual speckle noise level in the image (e.g. Marois et al., 2008; Wahhaj et al., 2013). Contrary to the classical SNRt maps, these SNRt,opp maps now suggest that parts of the forward-scattered edge of the outer disc are significant (Figure 10c and d). Nonetheless, this procedure still suffers from the small number of independent apertures at short radial separation.

An alternative way to assess the significance of signals in HCI post-processed images is to use standardized trajectory intensity mean (STIM) maps (Pairet et al., 2018a). STIM maps leverage the temporal variation of pixel intensities resulting from the trajectory of speckles in the derotated cube of residual images (i.e. the cube of images obtained after PCA modeling and subtraction, and subsequent alignment with North up and East left). These 2D detection maps are defined as $\mu(x_{i,j})/\sigma(x_{i,j})$, where $\mu(x_{i,j})$ and $\sigma(x_{i,j})$ are the mean and standard deviation of trajectory $x_{i,j}$ (i.e. the transversal slice at location ($i$,$j$)) throughout the derotated cube of residual images, respectively. Pixels containing authentic circumstellar signals (disc or planet) present a significantly higher value of $\mu(x_{i,j})/\sigma(x_{i,j})$ than pixels with no circumstellar flux contribution, due to (1) the relatively more constant (and higher) intensity of pixels containing circumstellar signal, and (2) the relatively higher standard deviation for pixels devoid of circumstellar signal due to the trajectory of speckles crossing these pixels in the derotated cube. By construction, STIM maps address the two main drawbacks of applying SNRt maps to our PDS 70 images: the presence of bright circumstellar signals spanning a range of radii and small number of independent apertures at short separation. The original definition of STIM maps was made in the context of 3D ADI datacubes (Pairet et al., 2018a). However, it is conceptually similar to apply it to SDI data; the trajectories are only radial instead of azimuthal. Our STIM maps for PCA-SADI and PCA-ASDI are thus computed on the (3D) cube of residual images obtained after the second PCA (ADI or SDI, resp.), for which the only change compared to the original STIM map definition is that it benefitted from a better whitening of the noise than using only a single PCA-ADI or PCA-SDI. They are shown in Figures 1c and d, and 10e and f.

In order to interpret the STIM maps and identify significant features, it is necessary to compute a STIM map using opposite parallactic angles (hereafter referred to as inverse STIM map). As explained for the SNRt,opp map, the inverse STIM map is expected to be a good proxy of the residual speckle noise in the (correctly derotated) STIM map. Figure 11 shows the cumulative fraction of pixel values in the inverse STIM maps for PCA-SADI ($n_{\rm pc}^{\rm SDI}=1$; $n_{\rm pc}^{\rm ADI}=1$) and PCA-ASDI ($n_{\rm pc}^{\rm ADI}=7$; $n_{\rm pc}^{\rm SDI}=1$). Pairet et al. (2018a) suggested to use this kind of plot to identify the threshold pixel value above which signal is very unlikely to stem from the random combination of residual speckle in the (correctly derotated) STIM maps. In view of the cumulative fraction distributions, we consider threshold values of 1.1 and 1.8 for the PCA-SADI ($n_{\rm pc}^{\rm SDI}=1$; $n_{\rm pc}^{\rm ADI}=1$) and PCA-ASDI ($n_{\rm pc}^{\rm ADI}=7$; $n_{\rm pc}^{\rm SDI}=1$) STIM maps, resp. (Figure 1c and d, resp.).

Appendix B Reliability of spectro-astrometry inferred by NEGFC

Wertz et al. (2017) proposed to estimate the astrometric errors associated to residual speckle noise by injecting a statistical number of fake companions at the same radial separation and contrast as the companion but at different PAs, and measuring the deviations between parameters used for injection and NEGFC-inferred parameters. We followed the same method and injected 42 fake companions in separate copies of our original cube where the protoplanet was subtracted using the best estimates of $r$, PA and contrast($\lambda$) inferred in Section 4.1.1. Given the bright forward scattered edge of the disc to the W and SW of the image, we only injected the fake companions between PA=-30° and PA=170°, at 10° intervals, at the same radial separation as the protoplanet. We carried out the NEGFC optimization for 8 values of $n_{\rm pc}^{\rm ADI}$ (from 3 to 10 included), for which the companion was visually conspicuous in the final image (Figure 2). We then found the optimal $n_{\rm pc}^{\rm ADI}$ that minimized the mean deviation between the true ($r$, PA) used for the injection and the NEGFC estimates for the 42 fake companions. Then, we considered the standard deviation of these 42 deviations as the astrometric uncertainty associated to residual speckle noise. We found uncertainties of 4.8 mas and 2.9°, for $r$ and PA respectively.

To ensure that the spectrum inferred for the companion is not distorted by the PCA-ASDIλ algorithm, we followed a similar procedure as above for the inferred contrast. Given the prior lack of knowledge of the optimal value of $n_{\rm pc}^{\rm ADI}$, we injected 42 fake companions at the median contrast as a function of wavelength inferred by NEGFC (PCA-ASDIλ) for the authentic companion for values of $n_{\rm pc}^{\rm ADI}$ ranging in $[3,10]$. Then, we applied NEGFC (PCA-ASDIλ) to estimate the spectrum of each of the fake companions. We considered either a median collapse or a variance-weighted average of residual images, and noticed a slight decrease of deviations between recovered and injected spectra using the variance-weighted average. For each binned spectral channel, we found the optimal $n_{\rm pc}^{\rm ADI}$ that minimized the median deviation between the true injected flux and the NEGFC estimates of the flux for the 42 fake companions. This optimal $n_{\rm pc}^{\rm ADI}$ was then used to infer the final contrast as a function of wavelength of the companion (Figure 8). The error bars in Figure 12 are the standard deviation of the deviations between the injected and estimated contrast($\lambda$) using the optimal $n_{\rm pc}^{\rm ADI}$ for each channel, and are also used in Figure 8. Figure 12 shows that the spectrum of fake companions inferred by NEGFC (PCA-ASDIλ) is consistent with the true injected spectrum, for all binned spectral channels.

Appendix C The ASDI+Br$\gamma$ and SADI+Br$\gamma$ algorithms

ASDI+Br$\gamma$ is similar to what is referred to as ASDI in Close et al. (2014) and SDI+ in Wagner et al. (2018), with the difference that we correct for algorithmic flux losses in each channel by measuring the throughput with injected fake companions. Our detailed procedure is as follows:

We first extracted a raw spectrum of the star using 1-FWHM apertures in the median $H$+$K$ spectral cube before ADI processing. No stellar emission line standing out from the continuum was found around the Br$\gamma$ line, but we found an absorption line slightly shortward (at $\sim 2.161\mu$m). 2. 2.

We defined our line channels as all spectral channels lying within 4 spectral FWHM from the center of the Br$\gamma$ line, where the assumed gaussian spectral FWHM of the line was varied between 7.5, 15 and 30 Å (i.e. we considered 6, 12 and 24 spectral channels). An identical number of continuum channels were considered directly longward from the last line channel (to avoid the stellar absorption line at shorter wavelength). 3. 3.

The continuum channels were rescaled both physically and in flux to match the size and flux of a weighted average of the line channels. The weights followed a gaussian centered on the Br$\gamma$ line, with the 3 tested values of line FWHM, and were normalized to 1. 4. 4.

We applied PCA-ADI in concentric annuli, on both the line and continuum channels, individually. 5. 5.

In each channel, we estimated the throughput of the algorithm as a function of radial separation. This was obtained at each radial separation with the injection of 7 fake companions at different azimuths, considering the throughput as the median throughput value out of those 7 flux ratio measurements (the algorithmic throughput is the ratio of the retrieved flux to the injected flux). 6. 6.

We divided the PCA-ADI frames by the throughput map computed in the previous step. 7. 7.

The Br$\gamma$ image, $I_{\textrm{Br}\gamma}$, was calculated as a weighted average of the PCA-ADI frames of the line channels (Figure 9a). The weights were the same as in step (iii). The continuum image, $I_{\textrm{cont}}$, was computed as the median of the PCA-ADI frames of the continuum channels (Figure 9b). 8. 8.

We obtained the ASDI+Br$\gamma$ image (Figure 9c) by subtracting the continuum image from the Br$\gamma$ image: $I_{\textrm{ASDI+Br}\gamma}=I_{\textrm{Br}\gamma}-I_{\textrm{cont}}$. 9. 9.

Steps (i) to (viii) were repeated varying the number of principal components ($n_{\rm pc}^{\rm ADI}$) between 1 and 10. A lower $n_{\rm pc}^{\rm ADI}$ led to more correlated residual speckle noise between spectral channels, but a larger $n_{\rm pc}^{\rm ADI}$ produced a better raw contrast in each PCA-ADI image. We determined the value of $n_{\rm pc}^{\rm ADI}$ that optimized sensitivity for Br$\gamma$ line emission at the location of the companion to be 7, as this value minimized the throughput-corrected noise at the radial separation of the companion in $I_{\textrm{ASDI+Br}\gamma}$.

By contrast, SADI+Br$\gamma$ is similar to what is referred to as ASDI in Rameau et al. (2015). Our detailed procedure was as follows:

We applied steps (i) to (iii) defined for the ASDI+Br$\gamma$ algorithm. 2. 2.

We calculated the pre-ADI Br$\gamma$ image as a weighted average of the line spectral channels. Similarly, the pre-ADI continuum image is computed as the median of the continuum spectral channels. 3. 3.

For each spectral cube of our observing sequence, the pre-ADI continuum image was subtracted from the pre-ADI Br$\gamma$ image. This led to a 3D (ADI) cube consisting of continuum subtracted Br$\gamma$ images. 4. 4.

PCA-ADI was applied to that cube in concentric annuli. We also applied PCA-ADI on the continuum and (non-continuum subtracted) Br$\gamma$ cubes. 5. 5.

We estimated the throughput of PCA-ADI as in step v of the ASDI+Br$\gamma$ algorithm. 6. 6.

We divided the PCA-ADI frames by the throughput map computed in the previous step. This led to a final SADI+Br$\gamma$ image (Figure 9f), and to both continuum and (non-continuum subtracted) Br$\gamma$ PCA-ADI images ($I_{\textrm{cont}}$ and $I_{\textrm{Br}\gamma}$, Figure 9d and e, resp.), all of which were throughput-corrected. 7. 7.

We repeated the process (steps i to vi) for different $n_{\rm pc}^{\rm ADI}$ until we found the value of $n_{\rm pc}^{\rm ADI}$ that minimizes the throughput-corrected noise. The optimal $n_{\rm pc}^{\rm ADI}$ for SADI+Br$\gamma$ was 2.