© The Authors 2023

1 Introduction

Direct imaging is a powerful technique for characterising the properties of exoplanets, such as their effective temperature, atmospheric composition, orbital motion, and top-of-atmosphere inhomogeneity (e.g. Skemer et al. 2014; Macintosh et al. 2015; Rajan et al. 2017; Samland et al. 2017; Wang et al. 2018; Biller et al. 2021). However, direct imaging of exoplanets is challenging due to the high contrast and small angular separation between the planet and its host star. Overcoming this challenge requires the development of world-class instrumentation, carefully designed observation strategies, and advanced post-processing algorithms. Dozens of exoplanets and sub-stellar companions have been imaged so far, such as the famous multi-planetary system, HR 8799 bcde (Marois et al. 2008, 2010), β Pictoris b (Lagrange et al. 2010), and 51 Eridani b (Macintosh et al. 2015). Surveys such as the SpHere INfrared Exoplanets (SHINE) project (Chauvin et al. 2017; Desidera et al. 2021; Langlois et al. 2021; Vigan et al. 2021), Gemini Planet Imager Exoplanet Survey (GPIES, Nielsen et al. 2013, 2019), B-star Exoplanet Abundance STudy (BEAST, Janson et al. 2019), Young Suns Exoplanet Survey (YSES, Bohn et al. 2020a,b, 2021), and Code for Orbital Parametrization of Astrometrically Inferred New Systems (COPAINS) pilot survey (Bonavita et al. 2022) are shaping our general understanding of the giant exoplanet and host star demographics.

Coronagraphs are important components of the instrumentation used in direct imaging. They suppress a significant amount of stellar flux while maximising the throughput of planet flux with a tradeoff between contrast, angular resolution, and working angles. The vector Apodizing Phase Plate (vAPP), the successor to the Apodizing Phase Plate (APP, Kenworthy et al. 2007), is a pupil plane coronagraph that uses liquid crystals and direct writing techniques to modify the phase of incoming light. This produces a dark zone in the coronagraphic point spread function (PSF) where high-contrast companions can be detected (Snik et al. 2012; Otten et al. 2014). Compared with focal plane coronagraphs, vAPPs are more stable to tip-tilt instability. Unlike the PSF of APP, the PSF of vAPP has a 360-degree coverage around the central star. Depending on its design, there are several types of vAPPs that produce different shapes of dark zones and PSFS (Doelman et al. 2021). Two common types of vAPPs currently installed on telescopes are the grating-vAPP (gvAPP180), which produces two 180° or D-shaped dark holes, and the double grating-vAPP (dgvAPP360), which produces a 360° circular dark hole.

Although vAPP coronagraphs have been developed for about ten years, a consensus on the optimal post-processing strategy for vAPP data has yet to be reached. In particular, the two complementary coronagraphic PSFs of gvAPP180 add an additional challenge to the data reduction. Several techniques before this work have been applied to the vAPP data. For example, Otten et al. (2017) exploited the symmetry of the two PSFs of gvAPP180. They reduced the first on-sky dataset of the gvAPP180 on MagAO/Clio2 by rotating, scaling, and subtracting one PSF from its complementary PSF. Sutlieff et al. (2021) applied three techniques to their HR 2562 dataset also obtained by the same gvAPP180. The first one is to join the two dark holes from the two complementary PSFs and implement classical angular differential imaging (cADI). The second technique, ADI + principal component analysis (PCA), also operates on the joined dark holes but uses PCA to model and subtract speckles before derotation. The last technique is flipped differential imaging (FDI + PCA) which uses the symmetry of the two complementary PSFs. It rotates the PSFs of one side by 180 degrees, implements PCA on them to build a reference PSF for the opposite PSFs, and then subtracts it. This technique is similar to the method used by Otten et al. (2017). They find that with these three methods, cADI produces the strongest detection of the companion. Otten et al. (2017) and Sutlieff et al. (2021) find that the symmetry between the two complementary PSFs produced by gvAPP180 is not as good as expected. Wagner et al. (2020) used dgvAPP360 on the Large Binocular Telescope (LBT) to observe the protoplanetary disk around PDS 201 and reduced the dataset with the package Karhunen–Loève Image Projection (KLIP, Soummer et al. 2012; Apai et al. 2016). Doelman et al. (2022) designed an ADI + PCA based technique to reduce the L band integral field spectroscopy (IFS) dataset of HR 8799 obtained from the dgvAPP360 mounted on LBT.

The above algorithms are all spatial algorithms. The recent emergence of a few time domain algorithms, such as wavelet-based temporal suppression for exoplanet detection (Bonse et al. 2018), exoplanet detection based on PAtch COvariance (PACO) (Flasseur et al. 2018), temporal reference analysis of planets (TRAP, Samland et al. 2021), half-sibling regression on exoplanet imaging (Gebhard et al. 2022), speckle space-time covariance in high-contrast imaging (Lewis et al. 2023), and PCA-Temporal (Long et al. 2023), has piqued interest in exploiting the temporal correlation of noise in high-contrast imaging. The TRAP algorithm is an inverse technique that reconstructs the planet signal and systematic noise simultaneously in a linear regression model. Unlike spatial algorithms that use local pixels from other images to build a reference spatial model, TRAP uses non-local but causally related pixels to build a reference temporal model. The planet signal of a single pixel over time is a positive transit curve modelled by simulating the planet trajectory due to the sky rotation in pupil-tracking observations. The systematic noise, including the stellar diffraction pattern, speckle, instrument scattered light, wind-driven halo, sky background, and other atmospheric effects, usually impacts a large region on the detector. TRAP reconstructs the systematic noise of the assumed planet transiting region, named the reduction region, from reference pixels that share the same underlying sources of systematic noise as the reduction region.

Because gvAPP180 produces two coronagraphic PSFs simultaneously, TRAP can select non-local pixels not only from the same PSF, which is the case for the conventional coronagraphic PSF, but also from the other PSF, which doubles the total number of available reference pixels. Thus TRAP can make the full utilisation of and explore different pixel combinations of the gvAPP180 coronagraphic data, instead of discarding partial information as several previously implemented spatial algorithms did, such as cADI. The previous algorithms based on the spatial symmetry of the two PSFs by Otten et al. (2017) and Sutlieff et al. (2021) are limited by the asymmetric noise in the two PSFs. TRAP may break this limitation if the noise is temporally correlated. Besides, applying a temporal algorithm that uses the causal information between pixels, such as TRAP, can reveal the temporal correlation of the systematic noise in high-contrast imaging, which in return can benefit the improvement of adaptive optics (AO) systems and coronagraphs such as vAPPs.

In this work, we apply the TRAP algorithm on two types of vAPP coronagraphic PSFs produced by dgvAPP360 and gvAPP180. Our main focus is to investigate the selection of reference pixels for vAPP coronagraphic data, especially for gvAPP180. We compare the photometric results obtained by TRAP with those obtained by previously implemented algorithms that demonstrated the best results on these datasets.

In Sect. 2 we provide a detailed description of the data employed in this study. We explain our methods of implementing TRAP on the dgvAPP360 and gvAPP180 coronagraphic PSFs, including the selection of reference pixels, in Sect. 3. We also reproduce the spatial algorithm with the best planetary detection results on the same dataset In Sect. 4 we present the results of different methods and compare their photometric measurements of the planets. We discuss the results in Sect. 5 and summarise this work in Sect. 6.

2 vAPP datasets

We used three on-sky vAPP coronagraphic datasets summarised in Table 1. They are the first few science observations obtained by the vAPP coronagraphs. Two of them contain real companions that can be used to compare the performance of different algorithms.

2.1 dgvAPP360 dataset of HR 8799

The first one is the IFS dataset of HR 8799 planetary system obtained by Doelman et al. (2022) from LBT with the L/M-band (3–5 μm) InfraRed Camera (LMIRcam, Skrutskie et al. 2010) on 18 September 2019. One dgvAPP360 was installed on LMIRcam and is also compatible with the Arizona Lenslets for Exoplanet Spectroscopy (ALES), an adaptive optics integral field spectrograph working at 2–5 μm (Skemer et al. 2015; Stone et al. 2022). The field of view (FoV) of ALES for this dataset is 2.″5 × 2.″5 and covers the three innermost planets of the HR 8799 system (from furthest to closest): c, d, and e. The plate scale of one spaxel is 35 mas. The planetary system was observed in the L-band from 2.74 to 4.29 μm with a spectral resolution of ~35, resulting in 100 wavelength channels. There are 1300 science frames in each channel with an integration time of 3934 milliseconds for each frame. There are also 99 background frames with the same integration time per frame as the science frames: the first 13 frames were observed after the first 100 science frames and all the last 86 frames were taken after the last science frame. The total observation time is 1.7 h with a field rotation angle of 85.°63. The seeing was 0.″8–1.″1 during the observation.

The coronagraph used in the HR 8799 dataset is the dgvAPP360 that has two liquid-crystal layers with specifically designed phase patterns that can suppress the stellar flux and form a circular dark zone surrounding the core of the central star (Doelman et al. 2020). The dark zone covers a 360° angle around the central star and the resulting PSF is similar to PSFs of the focal plane coronagraphs. Because of the double-grating design, the dgvAPP360 does not split the star light to the right and left circular polarisation states as the previously designed gvAPP180 does and produces only one PSF. The working angle of this dgvAPP360 is from 2.7 λ/D to 15 λ/D and is aimed to minimise the stellar diffraction halo close to the stellar PSF core. The inner working angle (IWA) and the outer working angle (OWA) of vAPPs are defined as the smallest separation and the largest separation where the throughput of the planet reaches half of the maximum throughput of the region where the contrast reaches the target contrast (Doelman et al. 2021). The raw contrast can reach 10⁻⁴ in the dark zone in L-band. A dgvAPP360 coronagraphic PSF example of the HR 8799 dataset is shown in Fig. 1.

Fig. 1 On-sky coronagraphic PSF from the dgvAPP360 mounted on LBT/ALES at 3.73 μm. The IWA and OWA are shown by the two dashed circles in red. The area between them is the dark hole of the PSF where the planet is expected to be located. The image size of ALES is smaller than the OWA. The colourbar shows the flux intensity after background subtraction.

2.2 gvAPP180 datasets of HR 2562 and Altair

Two datasets come from the 6.5 m Magellan Clay telescope with the adaptive optics (MagAO) system at Las Campanas Observatory, Chile (Close et al. 2010; Morzinski et al. 2014). They were collected by the same instrument at 3.94 μm with a band width of 90 nm: the Clio2 Narrow near-IR camera (Sivanandam et al. 2006; Morzinski et al. 2015) and the gvAPP180 coronagraph (Otten et al. 2017). The detector has a size of 1024 × 512 pixels with a pixel scale of 15.85 mas pixel⁻¹. One dataset was obtained on 07 February 2017 (Sutlieff et al. 2021) targeting at the F5V star HR 2562 which hosts a brown dwarf companion, HR 2562 B, at a separation of 643.8 ± 3.2mas (Maire et al. 2018). This dataset used an ABBA nodding pattern: the stellar PSF was alternately imaged on the left (A) and right (B) half of the detector. There are 354 cubes with 10 frames in each cube, 3540 frames in total (A: 1770 and B: 1770), and each frame has an integration time of 4 s with a total on-target integration time of 3.93 h. The derotator was on for the first 45% of the frames and switched off for the remaining part of the observation, producing a total field rotation angle of 42.°29. The seeing was 0.″5–1.″3 during the observation. The target of the other dataset is a much brighter A-type star, Altair, one of the first on-sky datasets obtained using this gvAPP180. The science PSFs were imaged only on the left half of the detector without any nodding pattern. 300 science data cubes were taken, each with 20 subframes. The integration time of each subframe is 1 s and the total on-target observation time is 6000 s. The telescope derotator was off during observation, producing a total field rotation angle of 42.°05. There are also 20 sky frames with the same integration time as the science frame, 1 s per frame. The star PSFs of the science frames were saturated and another 7 science data cubes (7 × 20 = 140 subframes) and 4 sky frames were taken as the calibration data. The integration time per calibration frame was 280milliseconds. The seeing was 0.″4–0.″7 during the observation.

This gvAPP180 coronagraph produces two coronagraphic PSFs by splitting the two circular polarisation modes of the star light (Otten et al. 2017). There is also a leakage term between the two PSFs, an intact stellar PSF without any dark holes. The intensity of the leakage term is faint and we did not make use of it in this work. All objects in the field of view receive the same phase adjustment by the gvAPP180 coronagraph. Hence the planetary PSF has the same shape as the stellar PSF in each coronagraphic PSF. If there is a planet located in the dark hole of one PSF, it also exists at the same location in the bright side of the other PSF. The working angle of this vAPP is between 2 and 7 λ/D and works from 2 to 5 μm. A gvAPP180 coronagraphic PSF example of MagAO/Clio2 is shown in Fig. 2. Each PSF has a half dark side where the starlight is significantly suppressed and a half starlight-unsuppressed bright side. Thus we have four different regions from the two opposite PSFs (upper PSF and lower PSF): upper dark, upper bright, lower dark, and lower bright.

2.3 Object property

HR 8799 hosts four giant, wide-orbit planet-mass objects discovered by direct imaging (Marois et al. 2008, 2010). The FoV of our dataset only covers the three inner planets, HR 8799 cde. HR 8799 d and e have a spectral type of L6-L8, and HR 8799 c is likely to have a later spectral type (Zurlo et al. 2016; Bonnefoy et al. 2016). HR 2562 B is a substellar companion orbiting HR 2562 and has a spectral type of L7±3 (Konopacky et al. 2016). Altair is a bright A type star without any known companions. Table 2 lists the properties of the central stars and their companions from the literature.

Fig. 2 On-sky coronagraphic PSF from the gvAPP180 mounted on MagAO/Clio2 at 3.94 μm. One PSF (upper PSF) is above another PSF (lower PSF) with a leakage term in the middle. The red circles in dashed lines show the IWA and OWA. The dark holes of the two PSFs complement the FoV of each other. The four different regions of the two PSFs we defined in this work are marked with white text. The colourbar shows the flux intensity after background subtraction.

3 Post-processing algorithms

The main goal of this work is to implement and test TRAP on vAPP coronagraphic data. The unconventional coronagraphic PSFs of vAPP have more flexibility in choosing non-local training data for the systematics model compared to conventional ones. This approach provides a unique opportunity to investigate the noise spatial and temporal correlations in vAPP coronagraphic PSFs and optimise the planet detection using a combination of training data. Additionally, this study may shed light on the underlying systematic noise and contribute to the development of future unconventional PSF designs. We extend the TRAP pipeline from Samland et al. (2021) to support the dgvAPP360 and gvAPP180 coronagraphic data and increase options in reference pixel selection. We also reproduce the previously implemented algorithms on vAPP coronagraphic data. Algorithms implemented on all datasets are described below.

3.1 Algorithms implemented on the dgvAPP360/HR 8799 dataset

The coronagraphic PSF of dgvAPP360 is a single symmetric PSF with a circular dark zone. We can directly apply ADI-based algorithms to this dataset as a regular pupil-tracking dataset without additional processing of the PSFs. Two algorithms are implemented on this IFS dataset and each wavelength channel is reduced independently.

Custom speckle subtraction algorithm (Custom): The sparsely sampled background frames, PSF drift (~2 spaxels or 70 mas), and global patterns probably introduced by the ALES pipeline make it very challenging to apply standard algorithms such as ADI or PCA to this dataset. Therefore, Doelman et al. (2022) design a custom spatial-domain method based on ADI and PCA to reduce the data. We follow their description to implement it on this dataset. Every four science frames were binned. Ten principal components of the 99 background frames were fitted to and subtracted from the science frames. The residuals were combined as a stellar PSF model. Then the stellar PSF model was scaled and subtracted from the original science frames. Afterwards, the star-subtracted frames were divided into eight groups in chronological order. For frames in each group, the background was removed by fitting ten principal components of frames that are at least 15-min away from this group (~0.85 λ/D for HR 8799 e) to the frame and subtracted from it. Then a low-order polynomial fitting along every row and column was used to remove the global structure, which varies between frames. The residual frames were aligned to north and then median combined. To remove the residual low-frequency noise, the final image was convolved with a Gaussian 2D kernel with a standard deviation of five pixels and subtracted by it.

Fig. 3 Planetary trajectories and reference pixel selection of the dgvAPP360 dataset. Left panel: planetary trajectories of HR 8799 cde. The image drift during the observations is also included in the trajectory. Right panel: reference pixels for a supposed position (cyan asterisk) used in TRAP, as shown by the white pixels. The known planet signals are masked.

TRAP

TRAP is a forward model that simulates planet positions in the original dataset and reconstructs all shared systematics, including the background, stellar speckles and other noise together with the planet signal in science frames using a causal time series linear regression model. A set of basis vectors are constructed by the principal analysis of the light curves of reference pixels that are impacted by the same noise sources as the reduction region. Therefore, the problems of this dataset are solvable in TRAP. Furthermore, because TRAP is a temporal algorithm, it works better at shorter integration time, for instance, a few seconds when the speckle temporal correlations are not averaged out (Samland et al. 2021). Therefore, the 1300 science frames were directly put into TRAP without any binning operation and pre-alignment. We used the same PSF model as the custom algorithm. For the reduction region, a PSF model with a radius of 1.1 λ/D was used to build the planet trajectory and the transit light curve of each pixel. For the reference pixel selection, we adopted a similar selection rule as in Samland et al. (2021): an annulus at the same separation of the reduction region, pixels surrounding the reduction region, and symmetrical pixels of the reduction region on the opposite side of the star. Pixels affected by the real planet signal were excluded from reference pixels, which is the known companion mask. The annulus width of the reference pixels was set to be 20 pixels to ensure the number of reference pixels was enough to reconstruct the systematic noise. Figure 3 shows the actual planetary trajectories due to field rotation and an example of the reference pixel selection of a supposed planet position. We also excluded reference pixels within the IWA when reducing regions outside the IWA, which resulted in a higher S/N.

3.2 Algorithms implemented on the gvAPP180/HR 2562 dataset

The D-shaped PSFs of gvAPP180 require specially designed data reduction techniques different from those for conventional circular coronagraphic PSFs. Sutlieff et al. (2021) apply three techniques to the gvAPP180/HR 2562 dataset and find that cADI produces the strongest detection of the companion among the three techniques. We follow their method description and find similar results. Hence, we choose cADI as the conventional speckle subtraction algorithm to compare with TRAP.

The star drifted linearly towards the lower left of the detector when the derotator was on, about 30 pixels (475 mas) on the y-axis and 35 pixels (555 mas) on the x-axis. Centring on the PSF core, we cropped the upper and lower PSFs separately with a frame size of 170 × 170 pixels (2.″7 × 2.″7), covering the whole working region of the gvAPP180. Then we did the linear, dark frame and flat field correction. We subtracted frames of one nod from the other to remove the thermal background. Furthermore, we binned every five frames in time to increase S/N and save computational time. The final total number of frames is 708 from the two nods.

Classical ADI:

The two dark holes from the lower and upper PSFs were joined together without overlap, and then regions inside the IWA and beyond the OWA were masked. As a reference PSF, the median of all frames was subtracted. Then the residual frames were derotated to north. We took the median of the derotated frames as the final residual map of cADI.

TRAP

The large drift of PSFs in the HR 2562 dataset makes the background contain not only temporal variations but also spatial variations due to the drift, for instance, the change of patterns spatially fixed in the detector. Therefore, the background could not be well reconstructed by TRAP in the aligned images. So we put the background subtracted and aligned images into TRAP. The aligned but background unsubtracted frames were used as the variance dataset for uncertainty calculation. Additionally, we implemented TRAP on cube A and cube B separately as they belong to different nod positions and hence very different regions in the detector, sharing few common temporal variations. By visual inspection, cube B was more heavily influenced by the detector and instrument effects, including a dark reflection ghost on the bright side of the upper PSF and a bright spike of scattered light across the lower PSF. The median of one cube was taken as the PSF model, and thus we had an upper PSF model and a lower PSF model. Furthermore, the flux in the dark holes of the PSF models was set to zero so that pixels with negligible planet signal from this area are not included in the reduction region (simulated by the PSF model) and can be used as reference pixels instead. Because they are next to the reduction region, they share similar noise as the reduction region and can increase the accuracy of the systematics model. We found that a larger PSF model size had a better result for cube B due to the more diffused planetary flux in nod B. The final PSF model used to simulate the reduction region has a radius of 1.1 λ/D for cube A, while the radius is 2.1 λ/D for cube B.

The key question when applying TRAP for gvAPP180 is how to choose reference pixels for the reduction region. Based on the same selection rules for dgvAPP360, we added the ‘half annulus’ option for gvAPP180: splitting the annulus into dark and bright sides. We heuristically designed six choices to select reference pixels from the four different regions of the upper and lower PSFs:

• Upper dark: choosing reference pixels from the dark side of the same PSF;

• Upper bright: choosing reference pixels from the same PSF, but only from the bright side;

• Lower dark: choosing reference pixels from the dark side of the complementary PSF;

• Lower bright: choosing reference pixels from the bright side of the complementary PSF, but the same positions as the reduction region in the upper dark hole have to be masked since they also contain planet signals;

• Upper dark+bright (standard): choosing reference pixels from both dark and bright sides of the same PSF, the same as for a symmetric PSF (e.g. the dgvAPP360 PSF);

• Joined dark holes: choosing reference pixels from the two joined dark holes.

White pixels in Fig. 4 demonstrate the six reference pixel options when searching in the upper dark hole, where the companion is located. Pixels influenced by known companion signals in both upper and lower PSFs were excluded from the reference pixel selection. The selected pixel for demonstration in Fig. 4 is also the expected companion location and thus the reduction region overlaps the known companion mask. When reducing the position within the IWA, we used all pixels outside the reduced region within the IWA as reference pixels to ensure that the number of reference pixels is greater than the number of frames.

Fig. 4 Six reference pixel designs for an assumed planet position (cyan asterisk) in the upper dark hole of a gvAPP180 PSF, as shown by the white pixels. Upper dark: choosing reference pixels from the dark side of the same PSF; upper bright: choosing reference pixels from the bright side of the same PSF; lower dark: choosing reference pixels from the dark side of the complementary PSF; lower bright: choosing reference pixels from the bright side of the complementary PSF; upper dark+bright: choosing reference pixels from the dark and bright sides of the same PSF; joined dark holes: choosing reference pixels from the joined dark holes.

4 Results

4.1 Photometric measurement and spectra comparison of the dgvAPP360/HR 8799 dataset

The residual map of the custom algorithm and normalised S/N map of TRAP averaged in 3.52–4.10 μm are presented in Fig. 5. HR 8799 cde are clearly detected by the two algorithms. For the custom algorithm, negative planet signals were injected to extract the planet location and contrast to the star (Marois et al. 2010; Lagrange et al. 2010). The negative planet signals were injected at the assumed planet locations in each frame. Then the injected dataset was reduced by the custom algorithm. The injected planet location and contrast were optimised iteratively with the downhill simplex method (Nelder & Mead 1965) until the residual planet signals were minimised in the residual map (e.g. Stolker et al. 2019). For TRAP, we measured the photometry by fitting a 2D Gaussian model to the unnormalised S/N map, which is the contrast map divided by the uncertainty map, because fitting in the normalised S/N map may slightly bias the planet location. Instead of interpolating the contrast at this position (sub-pixel precision), we then reduced the best-fit position using TRAP to obtain the contrast and uncertainty. The uncertainty was then normalised by interpolating the normalisation profile calculated from the unnormalised S/N map.

Doelman et al. (2022) use bootstrapping to estimate error bars for the retrieved contrasts. Their contrast spectra have higher S/N and smaller scatter than our version of the custom algorithm due to the difference in detailed implementation. We adopt their results on behalf of the custom algorithm to compare with the results of TRAP. The stellar PSFs in the blue and red edges of the IFS wavelength channels are blurred, which might be caused by cross-talk between wavelength channels. We discarded these channels and reduced the remaining 81 channels covering 2.91–4.16 μm with TRAP. Doelman et al. (2022) only calculate the uncertainty for each channel from 3.36 to 4.14 μm where the dgvAPP360 has high transmission. For wavelength shorter than 3.36 μm, they combine multiple channels to obtain two representative points. Therefore, we compare the results only in the common spectral range. The planet-star contrast spectra extracted with the two methods are presented in Fig. 6. They have similar profiles and agree with each other within one sigma in most wavelength channels. The considerable scatter and low S/N between 3.3 and 3.5 μm are due to the absorption band of the dgvAPP360 coronagraph (Otten et al. 2017; Doelman et al. 2022). The contrasts (here refers to planetary brightness relative to the star) of the three planets increase towards longer wavelengths since the planets are brighter at longer wavelengths. The S/N curves of the three planets are consistent with the dgvAPP360 coronagraph transmission profile, peaking at 3.65–3.75 μm. In Table 3 we present the S/N and contrast averaged between 3.52 and 4.10 μm where the transmission of the dgvAPP360 is the highest. In this range, TRAP achieves slightly higher S/N and less noisy contrast spectra for HR 8799 c and d than the custom algorithm. For both methods, HR 8799 d is brighter than HR 8799 c and has a higher S/N, which is consistent with the literature contrast in the L′ band in Table 2. However, compared to the outer two planets, the contrast spectrum of HR 8799 e has stronger high-frequency noise by both methods. Though the custom algorithm achieves a higher S/N for HR 8799 e compared to TRAP on average, the deviation of its S/N between wavelength channels is also larger. This indicates some inconsistent noise at small angular separations between different wavelength channels.

4.2 The gvAPP180/HR 2562 dataset

We reduced the two MagAO/gvAPP datasets of HR 2562 and Altair using cADI and TRAP. Because HR 2562 has a known companion, we mainly used this dataset to explore the algorithms for the gvAPP180 datasets. The Altair dataset is used for further analysis of noise temporal correlations in Sect. 5.

Fig. 5 HR 8799 result maps of the custom algorithm and TRAP averaged between 3.52 and 4.10 μm. Left panel: the final residual map of the custom algorithm where the colourbar represents the residual level. Right panel: the normalised S/N map of TRAP where the colourbar represents the S/N level. The two maps should not be compared directly.

Fig. 6 Contrast spectra and S/N comparison of the custom algorithm and TRAP. The results of the custom algorithm are adopted from the published results in Doelman et al. (2022), covering 3.36 to 4.14 μm. Left panel: the extracted spectra of HR 8799 c, d and e by the two algorithms between 2.91 and 4.16 μm. Right panel: the contrast S/N. The light green profile in the right panel is the star flux normalised by its maximum value between 2.91 and 4.16 μm, which is an indicator of the transmission of the dgvAPP360 coronagraph. The grey shaded area shows the absorption band of the coronagraph, where the extracted planet signals are the worst. On average, TRAP achieves slightly higher S/N and smoother contrast spectra for HR 8799 c and d. The custom algorithm achieves higher S/N for HR 8799 e but with a larger deviation between wavelength channels than TRAP.

Fig. 7 S/N detection maps of the upper dark region of the gvAPP180/HR 2562 dataset reduced with the six reference pixel choices in TRAP. Cubes A and В were reduced separately. The location of HR 2562 В is highlighted by the white circles. North is up and east is left. The companion is clearly detected in both cubes A and В when selecting reference pixels only from the dark hole in the upper PSF, where the companion is located.

4.2.1 Reference pixel selection in TRAP for gvAPP180

At first, we introduce two definitions in the TRAP detection map resulting from the dark and bright sides of the gvAPP180 PSF: dark region and transition region. When considering the planet trajectory in the pupil tracking dataset of gvAPP180 projected into a single PSF, as in Fig. 2: it can either fall completely into the dark hole or cross both the dark and bright sides. If it falls completely on the bright side, this is equivalent to falling completely into the dark hole of the opposite PSF. We define the ‘dark region’ in the detection map as the area where all simulated planet trajectories fall only into the dark hole. As such there are upper and lower dark regions in the detection map that correspond to the dark holes of the upper and lower PSFs, respectively. The dark region is always smaller than the dark hole and should be distinguished from the ‘dark hole’. The rest of the detection map is referred to as the ‘transition region’, where all simulated planet trajectories fall into both the dark and bright sides of the PSF. The transition region calculated from the upper and lower PSFs covers the same area in the detection map. It spans twice the parallactic angle range. For the transition region, we selected reference pixels from the same PSF using the standard choice as used for the dgvAPP360 PSF in Sect. 3.1.

The location of HR 2562 В is in the middle of the upper dark region in the detection map. We reduced the upper dark region with the six reference pixel selections demonstrated in Fig. 4. The reference pixel annulus had a width of 15 pixels and the known companion HR 2652 В was masked by a mask with a radius of 11 pixels (~1.5 λ/D). The 12 S/N maps of the upper dark regions of cubes A and В are presented in Fig. 7. The S/N maps were normalised by a three-pixel-wide annulus in the upper dark region as a function of separation. We found that selecting reference pixels solely from the upper dark hole allowed for a clear detection of the companion in cubes A and B. The detection was more pronounced in cube A than in cube В due to the latter’s worse image quality. The bright patches visible on the left edge of the upper dark region in the detection map of cubes A and В are due to the starlight contamination from the bright side of the PSF. These features also appear in the residual map of cADI both in our reduction and that of Sutlieff et al. (2021, see their Fig. 2). But they do not invalidate the detection of HR 2562 В at the expected location based on Maire et al. (2018).

The above result suggests that we should only choose reference pixels from the same dark hole when reducing the dark region of HR 2652. Hence, we reduced the upper dark region with reference pixels from the upper dark hole and the lower dark region with reference pixels from the lower dark hole. To generate a synthetic detection map, we used the weighted average (Samland et al. 2021) to combine the contrast ω and its uncertainty σ from the two independent reductions of the transition region. For each pixel in the transition region: $$\omega =\frac{{\omega }_1{\sigma }_{{\omega }_1}^{-2}+{\omega }_2{\sigma }_{{\omega }_2}^{-2}}{{\sigma }_{{\omega }_1}^{-2}+{\sigma }_{{\omega }_2}^{-2}},\sigma =\sqrt{{\left({\sigma }_{{\omega }_1}^{-2}+{\sigma }_{{\omega }_2}^{-2}\right)}^{-1},}$$(1)

where σ₁ and σ₂ are the preliminary uncertainty before normalisation. The normalisation for the synthetic S/N map (ω/σ) was calculated by dividing the dark region by its robust standard deviation within the dark part of the annulus and dividing the transition region by its robust standard deviation within the transition part of the annulus as a function of separation.

Fig. 8 Normalised S/N detection maps of cADI and TRAP for HR 2562. The dark region is where all the planet trajectories fall completely into the dark holes of the gvAPP180 PSF, while the transition region is where the planet moves across the dark and bright sides in the PSF. North is up and east is left. The dark region of the upper PSF is in the northwest and the dark region of the lower PSF is in the southeast. The transition region of cADI is enlarged to also include the severe speckle contamination near the boundary of the dark and bright sides of the PSF. The expected location of HR 2562 В is highlighted by the white circle. The bright and dark patterns near the edge of the dark regions in the S/N map of TRAP are results of the imperfect boundary between the bright and dark sides of the gvAPP180 PSF.

4.2.2 TRAP and cADI comparison for gvAPP180

The final S/N map of TRAP for the HR 2562 dataset is the weighted combination of the two normalised synthetic S/N maps from cube A and cube В also with Eq. (1). To measure the astrometry of HR 2562 В in the TRAP reduction, we fitted a 2D Gaussian model to the weighted combined contrast map. The photometric measurement of the companion was obtained by reducing the optimised position from cube A and cube В using TRAP, and then by combining the contrasts and uncertainties weightedly. For cADI we used aperture photometry and cross-correlation to calculate the S/N map (Ruffio et al. 2018). To create a flux map, we cross-correlated the residual map of cADI with an aperture having a diameter of seven pixels (~FWHM). Each pixel in the flux map contains the total flux within that aperture in the residual map. The S/N map was obtained by normalising the flux map by its standard deviation within a seven-pixel-wide annulus as a function of separation. The transition region of cADI should also span twice the parallactic angle range, the same as the transition region of TRAP. However, there are still lots of residual noise patterns at the border between the transition region and dark region due to the speckle near the joined boundary of the upper and lower dark holes. Therefore we enlarged the transition region to include them. The transition region and dark region in cADI’s S/N map were also normalised separately by the part of annulus in their own area as a function of separation, the same as in TRAP. The companion contrast and position of cADI were also extracted by the negative planet injection method in the residual map. The S/N and error of the contrast were estimated by interpolating the S/N map at the optimised planet position.

The S/N detection maps of cADI and TRAP are compared in Fig. 8. TRAP has a stronger detection of HR 2562 В compared to cADI. The central excluded region of cADI is larger than the IWA due to cross-correlation used in the calculation. Table 4 lists the photometric measurement of HR 2562 В. Compared with cADI, TRAP increases the S/N by 54%. Both measured contrasts from our cADI and TRAP agree with the contrast in Sutlieff et al. (2021) within 1er, which is 3.05 ± 1.00 ×10⁻⁴. The contrast difference between our cADI and Sutlieff et al. (2021) is due to differences in implementation details, such as bad pixels correction, dark holes cropping and joining, planet position optimisation, and S/N calculation. But even if compared with the S/N measured by Sutlieff et al. (2021), 3.04, TRAP still increases the S/N by 43%.

5 Discussion

5.1 HR 8799 e in the dgvAPP360 dataset

In Sect. 4.1 we reported the spectral results of HR 8799 ede with the custom algorithm and TRAP. The innermost planet HR 8799 e has the lowest S/N and largest high-frequency noise in the contrast spectrum compared to the other two planets. This is due to its proximity to the IWA of the coronagraph, where the noise pattern changes significantly from the stellar PSF core to the dark hole. Moreover, only a small number of pixels are affected by the same underlying cause of systematic noise at this small separation, making it challenging to build either a temporal reference model or a spatial PSF model in this region. As a result, both the custom algorithm and TRAP have difficulty detecting planet e as effectively as the other two outer planets. This also explains why excluding reference pixels inside the inner working angle for TRAP when reducing the dark hole results in a higher S/N, since the noise patterns are quite different.

The parameters of the custom algorithm and TRAP were optimised for a specific wavelength channel, 3.81 μm, where the transmission of the dgvAPP360 is high. Optimising the custom algorithm for each channel to remove the wavelength dependent systematic noise may achieve smoother spectra. TRAP can also produce better results if the parameters (e.g. the principal component number used to reconstruct the systematics and reference pixel selection) are optimised for each wavelength channel. However, based on our experience with optimised parameters in 3.81 μm, the improvement of TRAP may not be obvious in this case because the results are insensitive to them once the parameters are close to the optimised parameters. It is also computationally expensive to optimise parameters for each wavelength channel for a total number of 81 channels in both methods. Since we wish to have a user-friendly and practical algorithm for vAPP coronagraphic data, optimising the parameters for each channel is beyond the scope of this work.

5.2 Systematic noise in gvAPP180 PSFs

5.2.1 Confirmation by injection experiments on the Altair dataset

In Sect. 4.2.1 we expected the dark holes of the upper and lower PSFs to share the most similar noise properties, since they are the only ones with starlight suppression. The joined dark holes were expected to be similar to a conventional coronagraphic PSF and thus share the most systematic causes. Therefore, we hypothesised that the joined dark holes would produce the best result. However, the best result was found choosing reference pixels exclusively from the dark hole where the companion is located, running counter to our initial expectations. To further investigate the noise behaviour in the four regions, we performed tests on the Altair dataset without any known companions. We injected a fake planet signal with a 5σ contrast at the middle separation of the upper dark hole and retrieved it using reference pixels from the four regions, upper dark hole, upper bright side, lower dark hole, and lower bright side, respectively. More details on the reduction process for this dataset can be found in Appendix A. We did 5 injections for each choice, and Table 5 presents the averaged retrieved results. Choosing reference pixels from the dark hole of the same PSF, the upper dark hole in this case, allows us to retrieve the signal with an uncertainty of within 1σ, as expected. It also has the smallest position deviation, which is within one pixel scale, indicating the uncertainty of the optimised planet position by TRAP. Choosing reference pixels from the bright side of the same PSF produces the second-highest S/N, ~3.5, while choosing reference pixels from the complementary PSF has the worst performance, regardless from the dark hole or the bright side. This suggests that while the gvAPP180 coronagraph significantly suppresses starlight in the dark hole, there are still underlying noise causes that are shared by the dark hole and bright side of the same PSF.

5.2.2 Noise correlation between the two complementary PSFs

The results of reference pixel experiments on both HR 2562 and Altair indicate that the noise varies significantly across the detector. In the two complementary PSFs, the noise caused by different sources is more pronounced than that resulting from the same sources. The same causes of systematic noise shared by the two complementary PSFs are the dynamic aberrations and (quasi-)static optical aberrations between the telescope pupil and the coronagraph, such as the imperfect AO correction for the disturbance of the Earth atmosphere and non-common-path aberrations between the wavefront sensor and the coronagraph (e.g. Martinez et al. 2013). After passing through the gvAPP coronagraph, the light path splits into two paths, which arrive at different two regions on the detector. This process introduces several factors that contribute to the systematic noise differences between the two complementary PSFs, including (quasi-)static optical aberrations after the coronagraph, instrument scattered lights, reflection ghosts, uncorrected response variance across the detector, and other artefacts (Morzinski et al. 2015; Otten et al. 2017; Long et al. 2018). Many of these factors are visible in the background-subtracted PSFs and affect different parts of the detector, exhibiting few homogeneous temporal variations. Therefore, TRAP cannot reconstruct the systematic noise for one PSF by choosing reference pixels from the complementary PSF.

The detection difference between cube A and cube B in Fig. 7 is another piece of evidence for the noise variance across the detector. While the best reference pixel choice is the same for both cubes, the extracted signal from cube B is generally lower than that from cube A because the signal in cube B is more dispersed. The distance between A and B nod positions is even larger than that between the two complementary PSFs, resulting in less correlation between the noise of cubes A and B. Visual inspection reveals that cube В is more affected by instrument scattered light and reflection ghosts on the detector, which explains why the extracted companion signal from cube В is lower than that from cube A. In fact, these artefacts also break the symmetry of the upper and lower PSFs, thereby degrading the performance of algorithms based on the symmetry of the gvAPP180 PSFs (Sutlieff et al. 2021).

Fig. 9 5σ contrast curves of the gvAPP180 datasets with TRAP. The contrast curve is calculated as five times the median of a three-pixel-wide annulus as a function of separation in the normalised uncertainty map. Left panel: the contrast curves of the upper dark region of cube A of the HR 2562 dataset reduced by choosing reference pixels exclusively from the dark hole (labelled as ‘dark’) or both dark and bright sides (labelled as ‘dark and bright’). The detection significance of HR 2562 В is also marked in the figure. Right panel: the contrast curves of the dark region of the Altair dataset reduced by choosing reference pixels only from the dark hole or both dark and bright sides. The contrast curve of the transition region is also compared in the right panel, which is not much worse than that of the dark region.

5.2.3 Noise correlation in the same PSF

We identified three sources of systematic noise that are still shared between the bright side and dark hole of the same PSF. Firstly, the wind-driven halo around the central star (Cantalloube et al. 2018; Madurowicz et al. 2019) keeps the noise at small angular separations still correlated on both dark and bright sides. Though we performed frame selections to remove frames strongly influenced by the stellar halo in the Altair dataset, there are still frames containing residual wind-driven halo extending to the dark hole. Secondly, two bright ghosts on the edge of the dark hole introduce noise from the bright side to the dark side. These wind-driven halo and ghosts are also discussed in Otten et al. (2017), and are common problems for the MagAO/gvAPP180. Thirdly, although the gvAPP180 can discard the majority of starlight in the dark hole, a small remnant of the starlight and the diffraction pattern still exist in the dark hole. The raw contrast of the MagAO/gvAPP180 is 10⁻⁵–10⁻⁴ in the dark hole (Otten et al. 2017), leaving the photon noise in the dark hole still higher than the noise in other empty regions of the detector. Some quasi-static wavefront aberrations induced by the instrument can also brighten the symmetric diffraction structure around the star and thus the corresponding part in the dark hole (Otten et al. 2017). Therefore, considerable part of the noise in the dark hole is still associated with the noise in the opposite bright side and thus shares similar temporal variations. These effects are particularly pronounced for bright primary stars, especially at small angular separations. So including pixels from the bright side into reference pixel selection can better model the systematic noise for bright stars. Whether one should include bright pixels for a faint star is a case dependent problem. As shown in the contrast curve of HR 2562 in Fig. 9, including pixels from the bright side can improve contrast at certain separations, but there is no monotonie relationship. While the contrast curve of Altair in Fig. 9 shows that including pixels from the bright side indeed improves the contrast at small separations, the advantage recedes at larger separations. When the primary star is faint, we suggest using the standard reference pixel setting (dark and bright pixels) for an initial blind search, and then only choosing dark reference pixels to test if it can achieve higher S/N for the companion if detected in the blind search. When the primary star is a bright star similar to Altair or the wind-driven halo is very strong during observation, using the standard reference pixel setting is recommended.

5.3 Effective detection area of gvAPP180 PSFs

The contrast curves of the transition region and dark region of Altair are also compared in the right panel of Fig. 9. The contrast is not much worse than that of the dark region, partly because the transition region is independently reduced twice and combined from the upper and lower PSFs. It suggests that TRAP is able to detect a bright companion if it is unfortunately located in the transition region of gvAPP180.

To compare the effective detection area between TRAP and cADI, we used the artificial planet injection method to compare the retrieved photometric measurement of the two algorithms. We took the 5σ contrast map obtained by TRAP as the baseline and injected the 5σ signals into the Altair dataset, covering 412–792 mas. Each time we injected the signal at one position, and reduced the data by cADI and calculated the S/N of that position. Figure 10 shows the measured S/N of all injected positions reduced by cADI. The retrieved S/N of cADI has a median of 2.5 in the dark region and is statistically less significant than 5σ. This is consistent with the results of HR 2562 that TRAP increases the S/N by 43% compared to cADI. The median of the retrieved S/N of cADI in the transition region is only 0.9, much less significant than 5σ. It shows that TRAP also outperforms cADI in the transition region and the relative improvement is higher than that in the dark region.

Although it is preferable to ensure that the planet is located in the dark region of the gvAPP180 coronagraph for ADI observations to achieve higher S/N, this is not always possible for blind planet searches when the planet location is unknown beforehand. One solution is to rotate the gvAPP180 coronagraph to make the dark region cover an angle of 360°, which requires multiple observations of the same target. By improving post-processing algorithms such as TRAP, it is possible to detect relatively bright companions such as brown dwarfs, with an angle coverage of 360° by a one-time pupil-tracking observation without rotating the gvAPP180 coronagraph.

Fig. 10 S/N retrieved with cADI of 5σ injection. The injected 5σ signal at each position comes from the TRAP 5σ map of the Altair dataset. Left panel: retrieved S/N map with cADI. The injected region is from 412 to 792 mas. The stronger detection at the inner part of the dark region is due to the higher false alarm probabilities caused by small sample statistics (Mawet et al. 2014). Only a small number of pixels are available for normalisation at small separations, and therefore the uncertainty is underestimated. Right panel: S/N histogram of the left panel. The blue line depicts the profile of the dark region with a median S/N of 2.5 and the orange line depicts of the profile of the transition region with a median S/N of 0.9.

5.4 TRAP on dgvAPP360 and gvAPP180 datasets

The dgvAPP360 and gvAPP180 coronagraphs produce different PSFs. We find that reference pixels from a symmetric annulus should be used for the dgvAPP360/HR 8799 dataset; for the gvAPP180/HR 2562 dataset only reference pixels in the upper dark annulus should be used; for the gvAPP180/Altair dataset reference pixels from symmetric annulus should be used. These reference pixel selections seem to differ between the dgvAPP360 and gvAPP180 datasets, but actually they follow the same selection strategy: choosing reference pixels from the same PSF as the supposed planet location. The dgvAPP360 only produces a single coronagrahic PSF, so reference pixels are chosen from this PSF; for the gvAPP180, the upper PSF should be used for the reduction of the upper dark hole and the lower PSF should be used for the reduction of the lower dark hole. A new gvAPP180 coronagraph has recently been commissioned on the Enhanced Resolution Imaging Spectrograph (ERIS) on the Very Large Telescope (VLT) and had its first on-sky results (Kenworthy et al. 2018b; Boehle et al. 2021; Kravchenko et al. 2022). If there are less artefacts, less static optical aberrations in the different path after the coronagraph and more homogeneous detector response in ERIS/NIX than MagAO/Clio2, reference pixels from the complementary PSF can also be used to reconstruct the systematic noise for the opposite PSF, which may increase the contrast at smaller separations where the available reference pixels are less. This is more promising for the future 39-m European Extremely Large Telescope (E-ELT), which will feature two gvAPP180 coronagraphs on the Mid-infrared ELT Imager and Spectrograph (METIS, Brandl et al. 2018) and Multi-AO Imaging Camera for Deep Observations (MICADO, Davies et al. 2016) instruments (Kenworthy et al. 2018a; Doelman et al. 2021).

5.5 Outlook

There is no perfect metric for comparing different algorithms as their underlying assumptions are different. Our comparison of TRAP, the custom speckle subtraction algorithm, and cADI is based on the photometric measurements of known companions to determine which algorithm can obtain the best contrast in the dark hole, rather than searching for potential companions in the field of view. Adopting metrics such as true positive rates, false positive rates, and F1-scores is helpful to further compare the detection maps of different algorithms (Cantalloube et al. 2020).

In Fig. 8 there are some dark and bright patterns near the border of the dark and transition regions in the S/N map of TRAP. These patterns result from the difficulty in modelling the systematic noise around the border of bright and dark sides. We used the half annulus in the dark hole as reference pixels, which is cut by a straight border between the dark and bright side, but the real border is more complex than a straight line. This ‘half annulus’ option is intended to best model the systematic noise in the middle of the dark hole where the planet is expected to be located. It is not an optimised way to search for companions near the transition region. How to optimally combine reference pixels from dark and bright sides to search for planets in gvAPP180 coronagraphic PSFs using TRAP is not covered in this work. One possible choice is to make the half annulus in Fig. 4 move with the reduction region and always surround the reduction region in the centre instead of fixing the half annulus on the dark side. This choice may improve the detection on the edge of dark holes. For the IFS data, our implementation of algorithms is independent between different channels. Adding spectral differential imaging (SDI) to ADI-based algorithms can be applied to improve the detection (Kiefer et al. 2021).

As few techniques have been tested on vAPP coronagraphic data, even after adding TRAP tested in this work, we recommend exploring more algorithms, such as locally optimised combination of images (LOCI, e.g. Lafrenière et al. 2007; Wahhaj et al. 2015), ANgular Differential OptiMal Exoplanet Detection Algorithm (ANDROMEDA, Cantalloube et al. 2015), forward model matched filter (FMMF, Ruffio et al. 2017) and PACO (Flasseur et al. 2018) to find the most robust algorithm for vAPP coronagraphic data. For gvAPP180 coronagraphic data, we suggest exploring more inverse problem techniques like TRAP. They may have better performance in the transition region than algorithms based on speckle subtraction because the inverse techniques track the planet signal in the data cube, estimate the maximum likelihood and some of them even do not require subtracting the reference PSF or derotating the science frames, which avoid introducing more artificial noise in post-processing procedures. With more and more installation of vAPP coronagraphs (Doelman et al. 2021), testing different algorithms on vAPP coronagraphic data is important to prepare the community for the future observations and data reduction using these coronagraphs. As they suppress the star and planet PSF in the same way and preserves the central stellar PSF in per exposure, these coronagraphs are intrinsically suitable for variability monitoring of exoplanets and stellar companions by direct imaging (Sutlieff et al. 2023).

6 Conclusions

By training a temporal model on non-local pixels, TRAP enables an investigation of noise-time correlation between the dual PSFs of the gvAPP180 coronagraph. We find that the systematic noise between the two PSFs of MagAO/gvAPP180 is quite different. Unlike the conventional PSF, the symmetric systematic noise patterns between the two complementary dark holes are not significant and are overwhelmed by noise resulting from different sources. On the contrary, the bright side still shares a non-negligible amount of noise from the same underlying sources with the dark hole of the same PSF, even though the starlight is significantly suppressed in the dark hole. Thus we should only choose reference pixels from the same PSF when reducing the dark hole. If there is less systematic noise of different sources between the two PSFs, the opposite PSF can also be included to the temporal model, which will improve detections especially at small separations. This is promising for new generation of instruments, such as VLT/ERIS/gvAPP180, and is also instructive on the designing of the gvAPP180 coronagraphs on E-ELT.

In this work, we implemented a temporal algorithm, TRAP, on three vAPP coronagraphic datasets of two types of vAPP coronagraphs: the dgvAPP360 on LBT/ALES and the gvAPP180 on MagAO/Clio2. We compared the performance of TRAP with previously implemented algorithms that have the best S/N in companion detection, an ADI + PCA based custom speckle subtraction algorithm and classical ADI. For the dgvAPP360 dataset, TRAP achieves consistent results as the custom algorithm from Doelman et al. (2022), with slightly higher S/N for the outer two planets HR 8799 c and d and a lower S/N for HR 8799 e. The contrast spectra of planet e obtained by both algorithms present stronger high-frequency noise compared to the other planets. For the gvAPP180 dataset, the contrasts retrieved by TRAP and cADI from Sutlieff et al. (2021) of HR 2562 B agree within one sigma, and TRAP improves the S/N by 43%. TRAP also has the potential to detect relatively bright companions in the transition region where the planet trajectory spans the dark and bright sides in the gvAPP180 PSF. It makes the effective detection area of TRAP cover an angle of 360° around the central star, while cADI can only have effective detections in the dark region. The recommended choice for reference pixels of TRAP for dgvAPP360 datasets is the same as the standard setting, mainly a symmetric annulus at the same separation of the planet. For gvAPP180 datasets, when reducing the D-shaped PSF, it should be the annulus in the same PSF, either the half annulus in the dark hole or the full annulus including both the dark and bright sides. We also encourage further testing of other post-processing algorithms on vAPP coronagraphic data, especially to achieve better detections near the IWA of the dgvAPP360 PSF and on the border between the dark hole and bright side of the gvAPP180 PSF.

Acknowledgements

We thank the referee for providing constructive comments that helped to improve the clarity of this paper. B.J.S. is fully supported by the Netherlands Research School for Astronomy (NOVA). J.L.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement No 805445. This paper includes data gathered with the 6.5 metre Magellan Telescopes located at Las Campanas Observatory, Chile, and the Large Binocular Telescope located in Arizona. The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporation partners are: The University of Arizona on behalf of the Arizona university system; Istituto Nazionale di Astrofísica, Italy; LBT Beteiligungsgesellschaft Germany, representing the Max-Planck Society, the Astrophysical Institute Potsdam, and Heidelberg University; The Ohio State University, and The Research Corporation, on behalf of The University of Notre Dame, University of Minnesota, and University of Virginia. We thank all LBTI team members for their efforts that enabled this work. We gratefully acknowledge the use of Native land for our observations. LBT observations were conducted on the stolen land of the Ndee/Nnee, Chiricahua, Mescalero and San Carlos Apache tribes. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), SciPy (Virtanen et al. 2020), Scikit-learn (Pedregosa et al. 2011), NumPy (Harris et al. 2020) and Matplotlib, a Python library for publication quality graphics (Hunter 2007).

Appendix A Reduction of the Altair dataset

The Altair dataset shares similar properties as the HR 2562 dataset: significant linear PSF drifts between frames and similar parallactic angle range. We used the same data reduction method of the HR 2562 dataset to pre-process and post-process the Altair dataset as described in Sect. 3.2. Because the central star Altair was saturated, we used cross-correlation to align the images. The star PSF model was created by scaling up the calibration PSF to match the integration time of the science frames. The frames were binned every 20 frames. A master sky background made by the median of the 20 sky frames was subtracted from the science frames.

To verify the result found for HR 2562 that choosing reference pixels from the dark hole of the same PSF is better than choosing reference pixels from other regions, we did five injections at a separation of 5 λ/D, the middle of the upper dark hole. Each time we injected a fake planet at a different angle and reduced the upper dark region using the four reference pixel choices: only from the upper dark hole, only from the upper bright side, only from the lower dark hole, and only from the lower bright side. The injected signal was the 5σ contrast of that position reduced by choosing reference pixels from the upper dark hole. We calculated the retrieved S/N, contrast deviation and position deviation to the injected values. We averaged the results of the five injections for each reference pixel choice and took the standard deviation as the 1σ error as presented in Table 5.

References

All Tables

All Figures

	Fig. 1 On-sky coronagraphic PSF from the dgvAPP360 mounted on LBT/ALES at 3.73 μm. The IWA and OWA are shown by the two dashed circles in red. The area between them is the dark hole of the PSF where the planet is expected to be located. The image size of ALES is smaller than the OWA. The colourbar shows the flux intensity after background subtraction.
In the text

Fig. 2 On-sky coronagraphic PSF from the gvAPP180 mounted on MagAO/Clio2 at 3.94 μm. One PSF (upper PSF) is above another PSF (lower PSF) with a leakage term in the middle. The red circles in dashed lines show the IWA and OWA. The dark holes of the two PSFs complement the FoV of each other. The four different regions of the two PSFs we defined in this work are marked with white text. The colourbar shows the flux intensity after background subtraction.

In the text

	Fig. 3 Planetary trajectories and reference pixel selection of the dgvAPP360 dataset. Left panel: planetary trajectories of HR 8799 cde. The image drift during the observations is also included in the trajectory. Right panel: reference pixels for a supposed position (cyan asterisk) used in TRAP, as shown by the white pixels. The known planet signals are masked.
In the text

Fig. 4 Six reference pixel designs for an assumed planet position (cyan asterisk) in the upper dark hole of a gvAPP180 PSF, as shown by the white pixels. Upper dark: choosing reference pixels from the dark side of the same PSF; upper bright: choosing reference pixels from the bright side of the same PSF; lower dark: choosing reference pixels from the dark side of the complementary PSF; lower bright: choosing reference pixels from the bright side of the complementary PSF; upper dark+bright: choosing reference pixels from the dark and bright sides of the same PSF; joined dark holes: choosing reference pixels from the joined dark holes.

In the text

	Fig. 5 HR 8799 result maps of the custom algorithm and TRAP averaged between 3.52 and 4.10 μm. Left panel: the final residual map of the custom algorithm where the colourbar represents the residual level. Right panel: the normalised S/N map of TRAP where the colourbar represents the S/N level. The two maps should not be compared directly.
In the text

Fig. 6 Contrast spectra and S/N comparison of the custom algorithm and TRAP. The results of the custom algorithm are adopted from the published results in Doelman et al. (2022), covering 3.36 to 4.14 μm. Left panel: the extracted spectra of HR 8799 c, d and e by the two algorithms between 2.91 and 4.16 μm. Right panel: the contrast S/N. The light green profile in the right panel is the star flux normalised by its maximum value between 2.91 and 4.16 μm, which is an indicator of the transmission of the dgvAPP360 coronagraph. The grey shaded area shows the absorption band of the coronagraph, where the extracted planet signals are the worst. On average, TRAP achieves slightly higher S/N and smoother contrast spectra for HR 8799 c and d. The custom algorithm achieves higher S/N for HR 8799 e but with a larger deviation between wavelength channels than TRAP.

In the text

Fig. 7 S/N detection maps of the upper dark region of the gvAPP180/HR 2562 dataset reduced with the six reference pixel choices in TRAP. Cubes A and В were reduced separately. The location of HR 2562 В is highlighted by the white circles. North is up and east is left. The companion is clearly detected in both cubes A and В when selecting reference pixels only from the dark hole in the upper PSF, where the companion is located.

In the text

Fig. 8 Normalised S/N detection maps of cADI and TRAP for HR 2562. The dark region is where all the planet trajectories fall completely into the dark holes of the gvAPP180 PSF, while the transition region is where the planet moves across the dark and bright sides in the PSF. North is up and east is left. The dark region of the upper PSF is in the northwest and the dark region of the lower PSF is in the southeast. The transition region of cADI is enlarged to also include the severe speckle contamination near the boundary of the dark and bright sides of the PSF. The expected location of HR 2562 В is highlighted by the white circle. The bright and dark patterns near the edge of the dark regions in the S/N map of TRAP are results of the imperfect boundary between the bright and dark sides of the gvAPP180 PSF.

In the text

Fig. 9 5σ contrast curves of the gvAPP180 datasets with TRAP. The contrast curve is calculated as five times the median of a three-pixel-wide annulus as a function of separation in the normalised uncertainty map. Left panel: the contrast curves of the upper dark region of cube A of the HR 2562 dataset reduced by choosing reference pixels exclusively from the dark hole (labelled as ‘dark’) or both dark and bright sides (labelled as ‘dark and bright’). The detection significance of HR 2562 В is also marked in the figure. Right panel: the contrast curves of the dark region of the Altair dataset reduced by choosing reference pixels only from the dark hole or both dark and bright sides. The contrast curve of the transition region is also compared in the right panel, which is not much worse than that of the dark region.

In the text

Fig. 10 S/N retrieved with cADI of 5σ injection. The injected 5σ signal at each position comes from the TRAP 5σ map of the Altair dataset. Left panel: retrieved S/N map with cADI. The injected region is from 412 to 792 mas. The stronger detection at the inner part of the dark region is due to the higher false alarm probabilities caused by small sample statistics (Mawet et al. 2014). Only a small number of pixels are available for normalisation at small separations, and therefore the uncertainty is underestimated. Right panel: S/N histogram of the left panel. The blue line depicts the profile of the dark region with a median S/N of 2.5 and the orange line depicts of the profile of the transition region with a median S/N of 0.9.

In the text