Free Access
Issue
A&A
Volume 663, July 2022
Article Number A53
Number of page(s) 25
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202040024
Published online 12 July 2022

© ESO 2022

1 Introduction

Thousands of exoplanets have been discovered in the past decades using a variety of observational techniques. Directly measuring the brightness of the planets using high-contrast imaging techniques, provides crucial information on their temperatures and chemical compositions, and therefore on their possible formation history. This is particularly achievable when using spectro-photometry in the near- and mid-infrared (e.g., Samland et al. 2017). In addition, it allows us to probe a unique parameter space, in particular at large orbital separations (typically >10 au for nearby stars). Several dedicated surveys aim to find young giant gaseous planets using the direct imaging technique (see Bowler 2016 for an overview of previous large imaging surveys). Such surveys take advantage of the fact that gaseous giant planets at an early age (<100 Myr) are bright at infrared wavelengths because of the energy released during the contraction of their atmospheric envelope (see, e.g., Marley et al. 2007). For this reason, many of the recent, present, and future direct imaging instruments focus on near- and midinfrared observations (e.g., SCExAO, Jovanovic et al. 2015; MagAO-X, Males et al. 2018; SPHERE Beuzit et al. 2019; GPI, Macintosh et al. 2014). In particular, the NaCo Imaging Survey for Planets around Young stars (ISPY, Launhardt et al. 2020) is a survey whose aim is to find new young (<100 Myr) substellar (≲72 MJup) companions at L′ filter (~3.8 μm) using the NaCo instrument (NAOS+CONICA, Rousset et al. 2003; Lenzen et al. 2003) at the Very Large Telescope (VLT), Paranal Observatory, Chile.

The main challenge in detecting and characterizing faint companions around bright stars is to reach very high contrast values in the immediate vicinity of the star. Using a coronagraph to block the stellar light has proven to be a reliable strategy, improving the contrast by several orders of magnitude (e.g., Chauvin 2018) as close as a few λ/D from the star (typically >100 mas with current planet imagers). In addition, as mentioned in Launhardt et al. (2020), the contrast between a planet and its host star is more favorable at mid-infrared wavelengths. However, the observations are still largely dominated by the stellar halo and the point spread function (PSF), obscuring the planet signal in the images, so several observational techniques were designed to tackle this challenge: angular differential imaging (ADI, Marois et al. 2006), reference differential imaging (RDI, Smith & Terrile 1984), and spectral differential imaging (SDI, Marois et al. 2000; Sparks & Ford 2002). Since ISPY uses only one filter and has no spectral information, it uses ADI as it allows the temporal subtraction of the PSF and the quasi-static speckles using the data of the objects themselves as reference. Since the sky rotates during the observing sequence, while the PSF and the speckles remain quasi-static, subtracting the median from the datacube essentially removes their contribution in the final image (classical ADI, Marois et al. 2006). The suppression of the stellar PSF is significantly improved by using advanced post-processing algorithms and, for the forward- modeling approaches, prior information on the noise statistics and the ADI planetary signal itself. Some of these algorithms are the locally optimized combination of images (LOCI, Lafrenière et al. 2007); Karhunen-Loéve image projection (KLIP, Soummer et al. 2012), principal component analysis (PCA, Amara & Quanz 2012), KLIP-forward modeling (KLIP-FM, Pueyo 2016); the maximum likelihood-based ANDROMEDA (Mugnier et al. 2009; Cantalloube et al. 2015), local low-rank, sparse, and Gaussian noise component decomposition (LLSG, Gomez Gonzalez et al. 2016); a dedicated PSF reference library construction (e.g., Xuan et al. 2018; Ruane et al. 2019; Bohn et al. 2019) used to perform reference differential imaging (RDI, Lafrenière et al. 2009; Soummer et al. 2011; Gerard & Marois 2016); exoplanet detection based on patch covariance (PACO, Flasseur et al. 2018, 2020); and more recently the temporal causal regression model (lightcurve) approach (TRAP, Samland et al. 2021).

Surveys with very large samples allow us to put statistically robust constraints on the occurrence of giant planets around young stars, and therefore provide critical inputs to inform and refine planet formation models (see, e.g., Vigan et al. 2021). As many bright stars (L < 6 mag) were included in the ISPY sample, it was possible to observe a large fraction of them (~70%) using the annular groove phase mask (AGPM) vector vortex coronagraph (Mawet et al. 2013; Absil et al. 2014). The position of the star is a crucial parameter in the post-processing of ADI observations since it is the reference center to align and combine all the images. However, as soon as the star is placed behind the AGPM, its exact position can no longer be measured accurately, and this represents a challenge for high-resolution imaging observations, especially for first-generation direct imaging instruments such as NaCo. Second-generation instruments tackle this problem by using the deformable mirror to produce diffractive attenuated copies of the stellar PSF, called satellite spots (SPHERE, Beuzit et al. 2019; SCExAO/CHARIS, Groff et al. 2015) or by using diffracted images of the star produced by a reference grid inside the optical path of the instrument (GPI, Macintosh et al. 2014), to determine the position of the star behind the coronagraph with subpixel precision (reaching 0.1–0.2 pixels in the case of SPHERE in H band1).

Ground-based observations at mid-infrared wavelengths are performed with a very short detector integration time (DIT) to avoid saturating the detector, because of the significant sky thermal emission. The typical on-source time per target for the ISPY observing strategy is between 2 and 4 h, leading to 10000 to 30000 frames per target. Because the observing conditions vary throughout the night, or the coronagraph is slightly drifting (see Sect. 3.1), the quality of these thousands of frames will significantly vary. It is therefore relevant to assess whether removing the worst frames while keeping the most homogeneous ones can improve the final reduction. Different frame selection techniques have been proposed to try and improve the contrast or increase the S/N of any point source detection at different wavelengths. For example, Stolker et al. (2019) used a selection based on the sky brightness, Keppler et al. (2018) used the PSF variation of the stars in the field of view, while Bohn et al. (2020) discarded images with poor AO correction. There are therefore several approaches to performing some sort of frame selection, but overall they are usually on the conservative side, and there has not been an in-depth study of the potential benefits of performing more aggressive frame selection on large datacubes.

In this paper we present a new technique that aims to accurately find the position of the star behind the coronagraph for the particular case of NaCo observations with the AGPM. A good determination of the position of the star is essential for astrometric studies and a better extraction of the flux of potential companions. We also present our approach to registering the quality of the science frames, and to studying the impact of performing frame selection on the final reduction.

The paper is structured as follows. Section 2 summarizes the observation strategy and data reduction. In Sect. 3 we explain the new centering technique and the determination of the position of the star behind the coronagraph. We also explain our approach to registering and selecting the science frames. In Sect. 4 we benchmark our centering technique using observations of a star with a known companion. We apply our pipeline to two known companions (βn Pictoris b and R CrA b), ad present in Sect. 5 the improvements in the S/N of these point sources. In Sect. 6 we summarize our findings and conclusions on the importance of centering and frame selection for ADI observations.

Table 1

Log of the observations.

2 Observations and Data Reduction

In this section we describe the observing strategy, the standard (cosmetic) reduction steps, and the post-processing used to obtain the final images. In addition, we describe the effective shape of the coronagraphic image in the corrected science data that is described later on for the centering algorithm and registration of frames.

2.1 Observational strategy

The stars HD 34282, RCrA, and HD 104237 were observed as part of the ISPY program (Launhardt et al. 2020), during the nights 2016 November 07, 2017 May 19, and 2017 May 16, respectively. The observations were performed under different weather conditions varying from good to bad seeing (ranging between 0.34″ and 1.54″), as summarized in Table 1. We used the L27 camera (pixel scale ~27.2 mas pix−1) with the L′ filter (λc = 3.8 μm, Δλ = 0.62 μm) and the AGPM coronagraph for all three stars. The observing strategy was already described in Launhardt et al. (2020), but as a brief summary, it consists in obtaining multiple datacubes of 100 images each with the star manually placed behind the AGPM (hereafter science sequence). At the end of each science sequence the field of view is shifted so that it does not include the star (hereafter sky sequence). The purpose of the sky sequence is to estimate the thermal background and subtract it from the science observations. Each science-sky sequence lasts between three and five minutes, and is repeated tens of times for each star. At the beginning and end of every observation, we also obtained a separate set of images with the star placed at different regions of the detector (far from the AGPM), with shorter exposure times to measure the stellar PSF for photometric calibration (hereafter photometry frames or images). Since the third quadrant (bottom left corner of the camera) contains persistent bad columns, it was not used for the photometry observations.

β Pictoris was observed on 2013 February 01 (observations not part of the ISPY program), with the AGPM, under rather poor conditions. The data were presented in Absil et al. (2013) and re-reduced in Stolker et al. (2019), using the same strategy as used within ISPY, but for those observations the DIT was set to 0.3 s and the number of exposures per science cube was set to 200. The flux calibration strategy was also slightly different: in ISPY the star is moved to the three good quadrants of the detector, but for those observations sky measurements were taken between the photometry frames. The observations were executed using a window size of 768 × 770 pixels, which, in combination with the DIT, lead to a loss of ~ 10% of the frames. The total number of science sequences, the exposure times and the weather conditions are summarized in Table 1 for all stars.

2.2 Cosmetic Corrections

To prepare the data for the analysis, it is necessary to perform cosmetic corrections. First, we corrected the images for the dark current subtracting the suitable master darks. Master darks were produced by calculating the median of the dark frames grouped by the different exposure times used in an observing sequence (photometry, science, sky, flat field). Every image was corrected by the corresponding master dark by subtracting it from the image. The second correction corresponds to normalizing the response of all pixels. To this end, we constructed a master flat field by median averaging all the sky flat images (flat images are usually taken once per observing run and correspond to sky flats), and dividing by the mean value of counts. The correction for the flat field is done simply by dividing every image by its respective master flat field. At the end, all science, sky, and photometry images were corrected by master dark and master flat field. We also created a bad pixel map using the master dark and performing σ-clipping to select the outlier pixels (with σ = 5). The bad pixels were then corrected for by using a bicubic interpolation for every image. In the end, the sky mean was computed for every sky datacube using all available images. Then, the science images were sky-subtracted using the sky observations that are closest in time. We note that performing the sky subtraction this way is possible because of the negligible motion of the circular aperture and of the AGPM between consecutive science-sky sequences (see Appendix A).

2.3 Effective Shape of the Coronagraphic Image

When we subtract the sky background from the science images, we are also removing the thermal contribution of the AGPM (see Absil et al. 2016) and dust and debris located in the optical path of the telescope and instrument. A bright torus becomes visible around the coronagraphically rejected central star image. This torus is the result of the finite bandwidth of the AO, both spatial and temporal. Wavefront aberrations outside this bandwidth are not corrected for and lead to off-axis modes that are not rejected by the AGPM. Although the AGPM efficiently cancels out the on-axis light from the star, its transmission rapidly increases for off-axis light (see Fig. B.1 in Launhardt et al. 2020). The shape of the torus strongly depends on the alignment between the star and the AGPM, the performance of the adaptive optics system, and in general, the weather conditions. Both the AGPM and the star can move during an observing sequence (see Sect. 3), with the motion of the AGPM being negligible for consecutive sciencesky frames (see Figs. A.1 and B.1, respectively). For this reason it is important to find the center position of the AGPM on the detector and to find the star position with respect to the AGPM center.

2.4 Principal Component Analysis

The PCA and the derotation processes were performed using the PynPoint2 package (Amara & Quanz 2012; Stolker et al. 2019). The PCA process consists of creating an orthonormal basis of images (or eigenimages) representative of the ensemble of frames. Ultimately, the goal is to build an accurate representation of the stellar PSF and quasi-static speckles, which will be subtracted from each individual frame. Ideally, due to the sky rotation during the pupil-stabilized observations providing angular diversity with respect to the speckle noise, any real astrophysical signal should survive the process. Once the estimated PSF and quasi-static speckles have been removed, each frame is derotated according to its parallactic angle and all the frames are median-combined to recover the signal in the field of view (for more details, see Amara & Quanz 2012 and Hunziker et al. 2018). For the method to work, there is a compromise between the total field rotation (effective rotational angle) and the minimization of the flux loss of the astrophysical signal. The final images are obtained using a specific number of principal components. Given that the computational time necessary for the entire process is proportional to the size of the datacube, we performed our tests using different image sizes, between 39 × 39 and 121 × 121 pixels. For a given star, several images are produced, with different numbers of principal components to find the maximum S/N of the point sources in the field of view (if any) and to study the impact of the frame selection process as a function of the number of principal components. Their number varies in the range 1–30 or 1–150 (depending on the size and total number of frames).

3 CenteR Algorithm

In this section we introduce a new algorithm, CenteR3, designed to try to improve the data reduction process for angular differential imaging observations, focusing on observations taken with a coronagraph (AGPM). Our approach is twofold: first, we aim to estimate more accurately where the star lies behind the AGPM; second, we perform frame registration to characterize the quality of the science frames in order to perform frame selection later on.

thumbnail Fig. 1

Example of a sky image showing the circular aperture in blue (observing date 2016 November 08. The top right inset shows a zoom-in on the central part, with the AGPM thermal emission in the center. The white elongated spots correspond to dust or dirt in the optical path of the instrument or on the camera, and are emitting at near-IR wavelengths.

3.1 Calibrating the AGPM Detector Location

One critical aspect when performing high-contrast imaging with ADI, is to accurately determine the position of the star behind the coronagraph. This is an especially challenging problem for first-generation instruments such as NaCo since Mawet et al. (2013) and Huby et al. (2015) showed that the location of the AGPM is not always the same as that registered on the detector, physically on top of a less stable and dynamic adaptive optics correction. We therefore present here a new method to determine the position of the star behind the AGPM. When the star is aligned behind the AGPM, the resulting shape is similar to a torus for sky-subtracted science images (see Sects. 2.3 and 3.2 for an example). Therefore, the most straightforward way to obtain the position of the star is to fit a positive and a negative 2D Gaussian profile, modeling the aforementioned torus profile. However, the problem is highly degenerate if all the parameters are free to vary (around 13, depending on the assumptions), leading to poor constraints on the position of the star (uncertainty on the order of 0.65 pixels when considering all the parameters as free). We present a novel approach that helps determine the position of the star by fixing some of the free parameters, hence alleviating the degeneracies of the modeling. In particular, we fix some parameters of the negative 2D Gaussian profile related to the AGPM.

Our approach consists in finding the position of the AGPM independently, to keep it fixed during the rest of the modeling of the torus later on. In the sky images the AGPM is bright (see Fig. 1) due to its thermal emission, and as noted by Absil et al. (2016), it directly corresponds to the vortex center, which is not contaminated by the contribution of the star. By fitting a positive 2D Gaussian profile, for example, it is possible to accurately determine its position. However, for the science frames, the star is behind the AGPM, making it more difficult to determine the position of the AGPM. In Fig. 1, the AGPM is clearly visible close to the center, along with the borders of a circle that correspond to the edges of the 15″ circular aperture of the coronagraphic mask. We can correlate the center of this circular aperture with the position of the AGPM in the sky images. It is then possible to find the position of the AGPM in the science images by measuring the center of this circular aperture, following these steps:

  1. Measure the AGPM position in the sky frames;

  2. Measure the center of the circular aperture in the sky frames;

  3. Relate both centers to obtain the transformation between the position of the AGPM with respect to the center of the circular aperture;

  4. Measure the center of the circular aperture in the science frames without correcting for the thermal sky emission;

  5. Use the transformation found in (3) to obtain the AGPM position in the science frames.

We use the sky images to estimate the relationship between the position of the AGPM with respect to the position of the center of the circular aperture. Using the collapsed frames helps increase the signal of both the aperture and AGPM, therefore minimizing the final uncertainty. This can only be achieved because of the lack of significant motion of both positions in one sequence (for ISPY, 100 frames per sequence). The AGPM position is fitted by collapsing the image along the X - and Y-axes to obtain a more robust solution for the position, increasing the signal of the AGPM and making the modeling faster (see Fig. 2). However, 2D fitting is also possible, and depending on the model used and the number of free parameters, uncertainties on the order of 0.07 to 0.15 pixels can be reached, while our approach yields uncertainties of 0.1–0.15 pixels. In addition, we find that the 1D and 2D models are both consistent. The solution using the 1D model is obtained by fitting two normal distributions using a nonlinear least-squares method, minimizing the following equation:

η=i=1M(Z1I1exp[ ξ1 ]I2exp[ ξ2 ]ZBkg)2.$\eta = \sum\limits_{i = 1}^M {{{\left( {{Z_1} - {I_1}{\rm{exp}}\left[ { - {\xi _1}} \right] - {I_2}\exp \left[ { - {\xi _2}} \right] - {Z_{{\rm{Bkg}}}}} \right)}^2}} .$(1)

The value Zi corresponds to the normalized sum of counts per column (or row), ξ1=(XiμAGPM)2/2σ12,ξ2=(XiμAGPM)2/2σ22${{{\xi _1} = {{\left( {{X_{\rm{i}}} - {{\rm{\mu }}_{{\rm{AGPM}}}}} \right)}^2}} \mathord{\left/ {\vphantom {{{\xi _1} = {{\left( {{X_{\rm{i}}} - {{\rm{\mu }}_{{\rm{AGPM}}}}} \right)}^2}} {2\sigma _1^2,}}} \right. \kern-\nulldelimiterspace} {2\sigma _1^2,}}{{{\xi _2} = {{\left( {{X_{\rm{i}}} - {{\rm{\mu }}_{{\rm{AGPM}}}}} \right)}^2}} \mathord{\left/ {\vphantom {{{\xi _2} = {{\left( {{X_{\rm{i}}} - {{\rm{\mu }}_{{\rm{AGPM}}}}} \right)}^2}} {2\sigma _2^2}}} \right. \kern-\nulldelimiterspace} {2\sigma _2^2}}$, and Xi the position of each column (or row); μAGPM corresponds to the position of the AGPM; σ1 and σ2 are the standard deviations of the two Gaussians, I1 and I2 the scaling factor for the two profiles, and ZBkg is the background level, all of which are free parameters. The variable M corresponds to the number of columns (or rows). The modeling uses two Gaussians to account for the thermal emission of the AGPM and its halo (see Fig. 2). The results show that the movement of the AGPM during a normal sky sequence of 100 images remains negligible.

To estimate the position of the circular aperture, we first identified the pixels that reside inside the aperture. To this end, we compute both the median absolute deviation (MAD, σimage) and the median using the entire image, and then for each column along the X -axis Xi (or row along the Y -axis, respectively), we count the total number of pixels Ni that are three times σimage above the median4. If we were to use a linear sampling along the X-axis, the wings of the distribution would be undersampled, which would then bias the fitting results. The sampling therefore follows a cosine distribution to more evenly distribute the number of points in polar coordinates. To find the coordinate XCA and radius RCA of the aperture, we then minimize the quantity D=i=1N(XiRCAcos(αi)XCA)2+(Ni/SCARCAsin(αi))2,$D = \sum\limits_{i = 1}^N {\sqrt {{{\left( {{X_i} - {R_{{\rm{CA}}}}\cos \left( {{\alpha _i}} \right) - {X_{{\rm{CA}}}}} \right)}^2} + {{\left( {{{{N_i}} \mathord{\left/ {\vphantom {{{N_i}} {{S_{{\rm{CA}}}}}}} \right. \kern-\nulldelimiterspace} {{S_{{\rm{CA}}}}}} - {R_{{\rm{CA}}}}\sin \left( {{\alpha _i}} \right)} \right)}^2},} } $(2)

where SCA is a scaling factor to account for slightly elongated shapes (i.e., the aperture is not perfectly circular) and αi = arctan[(Xi − XCA)SCA/Ni]. This transformation allows us to work in polar coordinates from the center of the aperture, and the fitting is performed for the X- and Y-axes independently. Figure 2 shows examples of fitting the AGPM position in a sky frame (left panel) and of the circular aperture (right panel).

Figure 3 shows the position of the center of the circular aperture for science and sky images (blue triangles and red circles, respectively) for the observations of HD 34282 and R CrA (left and right, respectively). For HD 34282, the observing sequence was interrupted, which explains the two separate blocks of points, indicating that the coronagraphic mask actually moves between observing sequences. The figure shows that the centers of the apertures for the science and the sky frames follow the same path, but also that the center of the circular aperture slightly shifts during the observing sequence.

For the sky frames only, Fig. 4 shows the position of the AGPM and the center of the circular aperture for the same three targets. It is important to note that the AGPM and the center of the circular aperture move following the same path in both cases. To estimate the relationship between the two positions, we make the simplest assumption (i.e., a constant value), which is evaluated for each observing sequence. We find no evidence of a more complex relationship between the two positions than this one (e.g., a time dependence; see Appendix B). These constant values are determined as

ΔX=XAGPMXCA¯±δ(XAGPMXCA),$\Delta X = \overline {{X_{{\rm{AGPM}}}} - {X_{{\rm{CA}}}}} \pm \delta \left( {{X_{{\rm{AGPM}}}} - {X_{{\rm{CA}}}}} \right),$(3) ΔY=YAGPMYCA¯±δ(YAGPMYCA),$\Delta Y = \overline {{Y_{{\rm{AGPM}}}} - {Y_{{\rm{CA}}}}} \pm \delta \left( {{Y_{{\rm{AGPM}}}} - {Y_{{\rm{CA}}}}} \right),$(4)

where Xagpm and Yagpm correspond to the measured position of the AGPM on the X- and Y-axes for every individual sky frame, XAGPMXCA¯$\overline {{X_{{\rm{AGPM}}}} - {X_{{\rm{CA}}}}}$ and XAGPMXCA¯$\overline {{X_{{\rm{AGPM}}}} - {X_{{\rm{CA}}}}}$ correspond to the mean value of the difference between the AGPM position with respect to the center of the circular aperture for the X- and Y-axes, and δ(XAGPMXCA) and δ(YAGPMYCA) are the standard deviations of these relative positions. For each science frame, knowing the position of the center of the circular aperture and using those two equations, it becomes possible to determine the position of the AGPM with a precision better than 1 pixel (on the order of 0.1 pixels). The accuracy of 0.1 pixels in the position of the AGPM contributes directly to the final uncertainty in the star position, and it is the same order of accuracy as second-generation instruments (e.g., SPHERE or GPI).

thumbnail Fig. 2

Examples of what the AGPM and the circular aperture look like in a sky image. Left: modeling of the thermal emission profile of the AGPM along the X- and Y-axes to obtain the coordinate center for a cleaned sky image. Right: modeling of the circular aperture along the X- and Y-axes to estimate the coordinates of its center in the same image. The colored area corresponds to a 95% confidence interval.

thumbnail Fig. 3

Position of the center of the circular aperture (with respect to the camera) during the observation sequence for sky and science images (red and blue, respectively), for HD 34282 (left), β Pictoris (middle), and R CrA (right); the axis scales are 1:5, 1:1, and 1:3, respectively. The grid represent pixels, and the colored crosses correspond to the uncertainty obtained in the fit.

thumbnail Fig. 4

Position of the center of the circular aperture (blue triangles) and AGPM (red circles) with respect to their median value, for sky frames only. Left: position for HD 34282. Middle: same, but for β Pictoris. Right: same, but for R CrA. The axis scales of the left, middle, and right panels are 1:8, 1:1, and 1:3, respectively. The colored crosses correspond to the uncertainty obtained on the fit.

3.2 Estimating the Stellar Location

With the position of the AGPM known, it becomes easier to determine the location of the star behind it. With our approach we can accurately determine the location of the AGPM if we know the center of the circular aperture. In Fig. A.1 we furthermore show that the latter does not significantly move during an observing sequence, meaning that we are able to know the location of the AGPM for all frames without assuming for instance that the AGPM is located at the local minimum in the sky-subtracted frames. The frames are cropped into 23 × 23 pixel images, centered on the AGPM position to minimize the computational time. Then, each cropped science frame is fitted with a model that has three different contributions to reproduce the torus: a positive 2D Gaussian distribution representing the imperfectly AO-corrected non-attenuated star image on the detector VStar, a negative 2D Gaussian function accounting for the cancellation of the on-axis star light produced by the AGPM Vagpm, and an estimate of the background level Z0. The last is estimated in an annulus outside of the AGPM between 5 and 11 pixels accounting for the halo generated around the AGPM, while the 2D (negative or positive) Gaussian profile takes the general following form: V=Aexp[ 12(Xμ)TR(θ)×Σ1×R(θ)T(Xμ) ](2π)2 Σ .$V = A{{\exp \left[ { - {1 \over 2}{{\left( {X - {\rm{\mu }}} \right)}^{\rm{T}}}R\left( \theta \right) \times {\Sigma ^{ - 1}} \times R{{\left( \theta \right)}^{\rm{T}}}\left( {X - {\rm{\mu }}} \right)} \right]} \over {\sqrt {{{\left( {2\pi } \right)}^2}\left\| \Sigma \right\|} }}.$(5)

Here A is the amplitude, µ is the center of the distribution (fixed to Xagpm and Yagpm for the AGPM), X are the pixel coordinates [x,y], ς is the covariance matrix (a 2D diagonal matrix with σx2$\sigma _{\rm{x}}^2$ and σy2$\sigma _{\rm{y}}^2$ along the diagonal, allowing for elongated profiles), and R(θ) is a rotation matrix that depends on the angle θ. For the distribution that represents the stellar contribution, the free parameters are AStar, µStar, ΣStar (both σx,Star and σy,Star), and θStar. For the AGPM, µAGPM is fixed to the positions found previously using the circular aperture, and the free parameters are the amplitude AAgpm and ςagpm (assuming that σx,agpm = σy,agpm, and therefore R(θ) and θAgpm is set to zero). One possible caveat is that for frames where the star is slightly off-center with respect to the AGPM, our assumption that the negative Gaussian profile is centered at XAgpm and Yagpm may not be entirely correct. Indeed, subtracting a negative profile in the wings of a positive profile may slightly shift the position of the minimum. Nonetheless, this remains a good first-order approximation, and our approach is further tested in Sect. 4.

When minimizing the χ2 to model the sky-subtracted science frames, they are first cropped to a size of 41 × 41 pixels, centered on the position of the AGPM. Furthermore, we impose that the variable σx/y,AGPM (i.e., the 2D Gaussian σ of the AGPM) cannot be greater than the value of σx,Star2+σy,Star2$\sqrt {\sigma _{{\rm{x,Star}}}^2 + \sigma _{{\rm{y,Star}}}^2}$. We also consider the following constraints: for both the star and the AGPM, µ and σ cannot be greater than the actual size of the cropped frame. For σagpm, its value is constrained from fitting the sky frames, and it can vary but cannot be greater than 2.5 pixels5. Considering the original set of 13 unconstrained free parameters, we have now 8 free parameters in our modeling. In the end, for each individual frame, we store the location of the center of the circular aperture, as well as the eight fitted free parameters, and their associated uncertainties. Examples of the negative and positive 2D Gaussian fitting are shown in Fig. 5, showing the sky-subtracted science frames, the models, and residuals. In the model images we also give the location of the AGPM, showing that it does not necessarily correspond to the location of the minimum flux. Overall, with our approach, we are able to accurately find the position of the star behind the AGPM in the sky-subtracted science images with typical uncertainties of 0.20 and 0.19 pixels on the X- and Y-axes (see Sect. 4).

thumbnail Fig. 5

Modeling example of sky-subtracted science images for HD 34282 and β Pictoris using both the negative and positive 2D Gaussian. For each sub-panel, the observations are shown on the left, while the best-fit model is shown on the middle, and the residuals in the right. The homogeneity (H) and the S/N of the torus (TorusS/N) is labeled for each frame. The color scale is the same for each subpanel. The white cross for each center panel marks the location of the AGPM.

3.3 Frame Registration and Selection

The motivation behind frame registration is to characterize every individual image to identify subsamples that worsen some aspect of the reduction in order to remove those frames in the datacube later on, performing frame selection. Depending on the objective, there are different ways of performing this registration and selection of frames that will affect the final noise distribution, the signal of the candidate, and/or its position. For example, when we are interested in obtaining the position and flux of a candidate, it would be ideal to be able to maximize the signal of the object and/or minimize the residual noise of speckles. For this, we can discard those frames whose speckles and stellar PSF/halo distributions are very different from the rest of the frames, for example. Another need could be trying to maximize the S/N of the candidate, where the idea is to minimize the noise at the same angular separation to be sure that the signal is not being considerably contaminated by speckles. To this end, it could be possible to use frames whose statistical properties of noise are very similar as a function of angular separation, discarding those whose distributions are considerably different. Or even when we only want to increase the signal of the object, eliminating frames where the signal of the object has been very diluted due to a poor AO correction or changing weather conditions. In practice, there is no unique way to do the frame registration and selection. The most common approaches are carried out either visually when the total number of images remains small (e.g., Samland et al. 2017) or using simple criteria such as the background level (see Stolker et al. 2019; Cugno et al. 2019), the AO performance and correction (Hagelberg et al. 2016), or by studying the photometric variability of the point sources in the field of view (Keppler et al. 2018).

To quantify the degree of similarity among frames it is necessary to characterize the torus shape described in previous sections. One way to do this is by measuring the distribution of speckles and inhomogeneities generated around the AGPM, produced by a misalignment between the star and the AGPM, the weather conditions, and/or a poor AO correction, for instance. Similarly, the ring-shaped contribution (torus with isotropic distribution) formed due to the star–AGPM combination, can give us information mostly about the centering. The best way to measure both is by using image reconstructors, such as the pseudo-Zernike moments. The pseudo-Zernike moments were defined in Teh & Chin (1988) as a tool to obtain information about images and their reconstructions. The pseudo-Zernike moments are defined from the (pseudo-)Zernike polynomials (Zernike 1934; Bhatia & Wolf 1954), and they comprise a complete set of orthogonal functions on the unit disk. Teh & Chin (1988) showed that pseudo-Zerinke moments are less sensitive to image noise than the conventional Zernike moments. They are routinely used in optics, atmospheric sciences, and atmospheric turbulence corrections using adaptive optics systems (see, e.g., Noll 1976; Ma et al. 2017; Fusco et al. 2006; or Sauvage et al. 2016).

One of the main benefits of using the pseudo-Zernike moment is its fast computational analysis when performing such decomposition, making it possible to analyze quickly over thousands of frames. We use the package IM (Rajwa et al. 2013) from R Core Team (2019) to obtain the pseudo-Zernike moments for a specific order n and repetition m on the sky-subtracted science frames, where n = 0, 1, 2, … ∞ and |m|n. The order n is related to the complexity of the reference image (pseudo-Zernike moment) and m to the azimuth distribution. We note that m = 0 is purely related to radial profiles, and therefore it is not representative of anisotropics in the reconstruction. With the coefficients obtained from the moment decomposition (hereafter ζ), it is possible to define two criteria: the homogeneity H and intensity (S/N) of the torus formed around the AGPM in the science frames. The homogeneity of an image is related to the distribution of the light around the coronagraph; an image is classified with a high H value when the brightness distribution is homogeneously distributed as a function of the azimuth. For each frame H is calculated combining the information of the second and third repetitions (denoted m = ±1 and m = ±2) of the pseudo-Zernike moments as H=1|n=15(|ζnm=±1|)|2+|n=15(|ζnm=±2|)|2.$H = 1 - \sqrt {|\sum\limits_{n = 1}^5 {\left( {|\zeta _n^{m = \pm 1}|} \right)} {|^2} + |\sum\limits_{n = 1}^5 {\left( {|\zeta _n^{m = \pm 2}|} \right)} {|^2}.} $(6)

Both repetitions are related to the azimuthal brightness distribution with respect to the AGPM; therefore, large values of both repetitions implies a rather inhomogeneous frame. The S/N of the intensity of the torus is calculated using the first repetition (m = 0) of the pseudo-Zernike moments as TorusS/N=n=15(|ζnm=0|2)σ(25<r<35),${\rm{Toru}}{{\rm{s}}_{{{\rm{S}} \mathord{\left/ {\vphantom {{\rm{S}} {\rm{N}}}} \right. \kern-\nulldelimiterspace} {\rm{N}}}}} = {{\sqrt {\sum\limits_{n = 1}^5 {\left({|\zeta _n^{m = 0}{|^2}} \right)}}} \over {\sigma \left({25 &lt; r &lt; 35} \right)}},$(7)

where r is the distance from the location of the AGPM in pixels. To estimate the S/N of the torus, its intensity is divided by the noise, which is the standard deviation measured in a ring between 25 and 35 pixels in radius. This region was chosen as it is free of speckles since we want a S/N measure that reflects only the intensity of the torus and the background noise of the image not affected by residual speckles. The numerator of the previous equation corresponds to the intensity of the torus, but preliminary tests suggest that it is not by itself an ideal parameter to classify the science frames. The first repetition m = 0 is purely related to the radial and isotropic distribution of rings, and consequently it is azimuthally homogeneous. Therefore, any inhomogeneity of the torus would not be reflected in TorusS/N. Figure 6 shows two examples of reconstructed images and their respective residuals using the pseudo-Zernike moments. In addition to the reconstruction using only m = 0, which corresponds to the torus intensity (without azimuthal components), the figure also shows the reconstruction using |m| = 1 and 2 corresponding to the homogeneity as well as the reconstruction using 1 ≤ |m| ≤ 12. Both values, H and TorusS/N, are also normalized between 0 and 1, for all the frames and for each object individually.

With the results from the 2D negative and positive Gaussian profile fitting and with the information from the pseudo-Zernike moments, it is possible to further inspect the parameter space to distinguish between the most homogeneous well-centered frames and the bad ones. Different combinations of parameters can be used to separate the good and bad frames. However, the parameter space strongly depends on the observing conditions. Therefore, we need to adapt the selection criteria for each dataset. We find that the most crucial parameters to identify the set of good frames are the following:

  1. The separation between the estimated locations of both the star and the AGPM;

  2. The S/N of the torus;

  3. The background level and its standard deviation;

  4. The goodness of fit of the torus (from the 2D positive and negative Gaussian fitting);

  5. The ratio of the estimated uncertainties to the corresponding best-fit values for each of the free parameters.

The separation between the star and the AGPM provides a direct estimate of the quality of the centering during the observations. The S/N of the torus can indirectly provide comparable information: we find that when the S/N of the torus is high, the star is usually well-centered and the surface brightness is homogeneously distributed around the AGPM (in comparison with a misaligned case under same conditions). However, the value of the S/N is also affected by the weather conditions (e.g., clouds), decreasing the overall signal of the torus even if it is symmetric around the AGPM. The background level and the background noise level around the AGPM give us information about the strength of the speckles produced (e.g., their magnitudes and radial distribution), and as a consequence the quality of the AO correction. A poor AO correction usually generates many speckles, thus decreasing the contrast close to the AGPM. The ratio of the uncertainty to the estimated parameter serves as a filter to identify the frames where the fitting procedure was not successful. Other parameters, such as H, provide similar information and we note that their incorporation in the analysis does not necessarily improve the frame selection process (see Appendix C). However, its use can be beneficial for the characterization of frames in instruments whose modeling of the torus can be very complicated or even when the torus is not clearly detected.

Our approach for the frame selection is shown in Fig. 7, for HD 34282, showing in the left panel the S/N of the torus as a function of the separation between the star and the AGPM positions along the X-axis. The right panel shows the separation between the star and the AGPM for both axes. For both panels, the points are color-coded according to the goodness of fit when modeling the torus with two 2D Gaussian profiles (i.e., the sum of the squared residuals within the 23 central pixels). The vertical cluster of points that appears on the left side of Fig. 7 can be explained as follows. When the star is well aligned with respect to the AGPM, the S/N of the torus increases as the starlight is more homogeneously distributed around the AGPM. We find that in these cases the fitting process yields better results (smaller residuals). When the separation between the star and the AGPM increases, the S/N of the torus decreases leading to higher values for the residuals. Therefore, the resulting distribution, as shown in Fig. 7, shows a relatively wide ring of points, with low S/N values, as well as a narrower cluster of points with higher S/N values. However, the torus is also affected by the weather conditions and the AO performance. Even if the star is well centered behind the AGPM, if the observing conditions degraded, the S/N of the torus will also decrease. The resulting distribution is therefore a mix between well or poorly centered frames, all of them affected by the observing conditions and AO performance.

In the left panel of Fig. 7 the horizontal plane formed around the normalized S/N value of 0.3 is most likely related to poor AO correction. When this happens during the observations the background level around the AGPM increases due to the stellar leakage, and more speckles appear close to the AGPM. The algorithm can then only poorly constrain the position of the star given the overall lower quality of the AO correction. It can even find the brightest speckle next to the AGPM as the best approximation for the star location (top right panel of Fig. 8). Therefore, the intensity of the torus is dominated by the stellar leakage background, which is much larger than for other frames as the stellar light is much more diluted and less affected by the AGPM obstruction. As a result, we observe this relatively constant level in the measured S/N when the AO correction is not working optimally. In these cases, we find that the speckles are randomly distributed in intensity and separation, and the fitting procedure finds separations up to several pixels between the location of the star and of the AGPM (instead of a fraction of a pixel).

The next step is to design the criteria that will be used for the frame selection. The goal is to select the most homogeneous frames of the sequence, the ones that share the most common features, which can then be removed using a principal component analysis (PCA), for instance. To this end, it is possible to reject the outliers in some of the distributions (e.g., star to AGPM separation, background level). This can be done by rejecting all the frames that are 5 × MAD away from the median. In addition, it is also possible to fine-tune the acceptance level for the S/N of the torus. This can be done by rejecting the frames that comply with being less than ϕ × MAD + τ, where τ is the median and MAD is the median absolute deviation for the S/N of the torus, while ϕ being a free parameter that can take values from –∞ to +∞, in principle. It is possible to also consider the ratio of the estimated uncertainty to the best-fit value for each of the parameters used when modeling the torus with the positive and negative Gaussians. For each frame, if the ratio is higher than ρ (being the second free parameter, with ρ between 0 and +∞), then the frame is not considered in the rest of the analysis.

Figure 8 shows an example of the frame selection for HD 34282 using ϕ = 0.5 and ρ = 0.5, keeping ~55% of the frames equivalent to 17 934 frames. The left panel shows the normalized TorusS/N as a function of the separation between the star and the AGPM location along the X-axis, with the red dots corresponding to the selected frames and the gray dots to the rejected ones. The histogram at the top of the left panel shows the distribution of the selected and rejected frames. The orange circles correspond to examples of selected frames, while the cyan triangles examples of rejected frames, and these examples are shown on the right side of the figure (bottom and top, respectively). The red crossed circles give the location of the AGPM, while the yellow stars are the estimated position of the star. As mentioned previously, the horizontal cloud of gray dots at 0.3 in TorusS/N corresponds to frames where the AO system was not working nominally.

For our frame selection method the most important parameter is ϕ as it directly relates to the percentage of frames that will later on be used in the PCA. Figure 9 shows an example of the number of selected frames as a function of ϕ for HD 34282. The parameter ρ helps us determine the fitted parameters that are poorly constrained and, according to Fig. 9, it has a rather marginal effect on the number of frames selected. For this reason, we prefer to use ρ = 0.5 as reference value.

thumbnail Fig. 6

Two examples of reconstructed images using the pseudo-Zernike moments. For each panel (made up of six subpanels): top left: observations; top center: reconstructed image using m up to 12; top right: residuals of the modeling; bottom left: reconstruction of the image using only m = 0, related to radial and uniform azimuthal distribution; bottom center: reconstruction of the image using |m| = 1 and 2, corresponding to the adopted inhomogeneous contribution; bottom right: reconstruction of the image using 1 ≤ |m| ≤ 12, related to all the azimuthal contributions. For each subpanel the color scale is the same.

thumbnail Fig. 7

Examples of the distribution and registration of frames as a function of different parameters. Left: S/N of the torus as a function of the difference between the positions of the star and of the AGPM along the X-axis. Right: relative position of the star with respect to the position of the AGPM. For both panels the star is HD 34282 and the color-coding corresponds to the goodness of fit. The crosses in both panels correspond to the typical uncertainties. The contours contain 93 and 68% of the data in the left panel, and 93, 68, and 50% of the data in the right panel.

thumbnail Fig. 8

Examples of selected or rejected frames after applying our selection criteria. Left: normalized TorusS/N as a function of the difference between the positions of the star and the AGPM along the X-axis. The color-coding corresponds to the frame selected using ϕ = 0.5 and ρ = 0.5 (red dots) and rejected frames (gray dots). The orange circles correspond to three examples of frames selected and the cyan triangles to examples of rejected frames. The small panel at the top of shows the histogram of both distributions in logarithmic scale. Right: three examples of rejected frames (top row) and selected frames (bottom row). The red crossed circles correspond to the position of the AGPM, while the yellow star gives the derived position of the star obtained from the negative and positive 2D Gaussian fitting. The color scale is the same for all six images.

4 On-sky Validation: Benchmarking the Centering Algorithm

In this section we present a test on the robustness of our method to recover the true position of the star behind the AGPM. Given that there is no direct and model-independent way to measure the true position of the star behind the AGPM for NaCo, the best way to validate our method is to use pair(s) of stars whose angular separations are known or measurable. Therefore, we analyze observations of a target with a bright nearby known companion. The central star is HD 104237 and its companion, HD 104237 B, which lies in the field of view of ~3.29 – 3.29′′ (or 121 × 121 pixels), is detected in all the individual science and photometric frames. Using the photometric frames, in which the central star is >1.6′′ from the AGPM center and unsaturated, we are able to measure the separation and position angle of the companion with respect to HD 104237. We can therefore use this information later on to estimate the location of the central star in the science frames, correcting the position angle for changes in the parallactic angle.

The angular separation and the position angle are measured as follows. First, we fitted 2D Gaussian profiles to obtain the positions of both point sources in the reduced photometric frames. This was done by cropping the frames to smaller windows and centering them on both components. To avoid being biased by the Airy rings or being limited by low number statistics, we chose different sizes when cropping the images (10 × 10, 15 × 15, and 20 × 20 pixels) and calculated the mean weighted by the respective uncertainties. Then we measured the angular separation and position angle, propagating the estimated uncertainties on the final values. For each photometric integration we used the median stacked image over the 400 frames, calculating the median parallactic angle and its range as uncertainty. The photometric sequence is performed three times at the beginning and three times at the end of the observations, leading to six measurements for the angular separation and position angle. However, since the first sequence was observed with a poor AO correction (the stellar flux was completely diluted in the image), we discarded these bad-quality frames. We obtained the following results: 50.69 ± 0.20 pixels and 254.633 ± 0.122° for the separation and position angle, respectively. With a pixel scale of 27.212 ± 0.096 mas pixel−1 and a true north of 0.394 ± 0.212° (see Launhardt et al. 2020), the projected separation measured corresponds to 1.379′′ ± 0.007 and position angle of 255.027 ± 0.245°.

HD 104237 (DX Cha) is an accreting Herbig Ae spectroscopic binary (A or AaAb, Böhm et al. 2004; Garcia et al. 2013) with a circumstellar disk. Four stars are surrounding HD 104237 A within 15′′, labeled with letters from B to E (Feigelson et al. 2003), forming a kind of mini-cluster (Grady et al. 2004). The five sources (A-E) are confirmed members of the ϵ Cha association (Murphy et al. 2013), while HD 104237 A is the only early-type star of the association (3–5 Myr, Feigelson et al. 2003). Grady et al. (2004) reported that HD 104237 B (“source B” in Feigelson et al. 2003, or “star-2” in Grady et al. 2004), is a M3-4 star with an angular separation of 1.365′′ ± 0.019 and position angle of 254.586° ± 0.350 with respect to HD 104237 A (observing date 2003 June 10). Our measurements (ISPY observations, 2017 May 16) of position angle and angular separation (Table 2) both agree well within the uncertainties with the values presented in Grady et al. (2004).

The next step is to estimate the position of the star behind the AGPM, in a model-independent way, using the position of the companion. For each individual science frame we measure the location of the companion by fitting a 2D Gaussian profile. With the parallactic angle, the position angle, and the angular separation we can infer the true location of the star behind the coronagraph without any assumption about the reference center. We note that during the observing sequence, there were some technical problems and the observations had to be stopped. As a consequence, there is a jump in the sequence, leading to two separate distributions. In addition, the observations were executed under variable weather conditions (see Table 1).

To estimate the robustness of our method in finding the position of the star using the 2D negative and positive Gaussian profiles, we performed two different tests. The first consists of directly comparing the position of the close companion with the previously inferred position of the star and AGPM (see Sect. 3). The purpose of this test is to estimate the impact of using either the star position or the center of the AGPM as the center of rotation. We derotate each frame to align the position of the close companion to the north using both its position angle and the parallactic angle of the frame. This is done in two separate ways: (i) using the AGPM as the rotation center and (ii) using the stellar position derived in Sect. 3.2 as the center of rotation. To center the 2D distribution around zero, we then compute the difference between the location of the companion and the location of the star (or AGPM) and subtract the angular separation of the companion measured in the photometric frames. Measuring the dispersion of the final density distribution allows us to test whether the center of rotation is the correct one. For this test we kept ~93% of the data as we only removed the frames for which we could not determine accurately the positions of the companion and of the star (uncertainties larger than 1 pixel). The results are shown in Fig. 10. The left panel shows the density contours when the derotation is performed centered on the modeled position of the star, while the right panel shows the results when centered on the position of the AGPM. The median absolute deviation correspond to 0.20 pixels for the X-and Y-axes using the star position as reference center. When using the AGPM as the reference center we measure dispersions of 0.22 and 0.27 pixels. If we assume that the chosen center is correct, we expect to obtain a distribution centered around zero with low dispersion. Since we derotated the companion to the north, dispersion on the X-axis is dominated by uncertainties on the parallactic and position angles, while the dispersion along the Y-axis is dominated by the uncertainty on the angular separation, as shown in the left panel of Fig. 10 where the position of the star was used as the center. However, when the position of the AGPM is used as a reference center, the distribution is elongated along the Y-axis, suggesting that the angular separation between the AGPM and the companion varies. We therefore conclude that using the AGPM as the rotation center can introduce additional biases in the reduction process, while using the position of the science target leads to a smaller dispersion of 0.20 pixels.

The second test consists in directly comparing the true position of the central star (inferred from the position of the companion) with respect to the position of the star obtained with our method, without performing any derotation, and therefore without having to assume a center. Figure 11 shows the relative position of the central star using the two methods. The MAD corresponds to 0.21 and 0.19 pixels on the X- and Y-axis, respectively. Compared with the first test, the results strongly suggest that our method agrees well with the true position of the star within a dispersion of 0.20 pixel.

From both tests we can therefore conclude that our approach to finding the position of the star behind the AGPM agrees with the true position of the star estimated using its companion (within an uncertainty of 0.20 pixel), and to the same order of accuracy as second-generation instruments (e.g., SPHERE or GPI). Huby et al. (2015, 2017) demonstrated that there are some nonlinear effects when the star is close to the center of the coro-nagraph installed for Keck/NIRC2, which we do not take into account with our approach. These effects can introduce a slight shift or increase the dispersion in our measurements on the position of the star. However, we do not observe this effect for NaCo (see Fig. 11), or at least it is not measurable considering our uncertainties.

thumbnail Fig. 9

Percentage of selected frames as a function of the free parameter ϕ for HD 34282. The different lines show the behavior for three different values of ρ.

Table 2

Astrometric measurements of HD 104237 B.

thumbnail Fig. 10

Density contours of the position of the close companion after derotating each frame by the parallactic angle and the position angle of the companion, and after subtracting its angular separation. The derotation is either done at the position of the star inferred from our modeling (left) or at the position of the AGPM (right). The black circle corresponds to the uncertainties on the angular separation and the position angle of the companion at the peak of the distribution. The grid shown in each panel represent pixels, and the colored crosses correspond to the typical uncertainty. The contour lines in each panel enclose the areas containing 99.7, 95.5, 68.5, 50, 30, and 10% of the sample with respect to the maximum of the distribution.

thumbnail Fig. 11

Comparison of the star positions derived from the close companion and those obtained by our fitting approach (negative and positive 2D Gaussian fitting labeled NP2DG in the figure). The small top and side panels show the kernel density estimators (KDE) for the X- and Y-axes, respectively. The grid represent pixels, and the colored cross corresponds to the typical uncertainty. The contour lines enclose the areas containing 95.5, 68.5, 50, 30, 10, and 1% of the sample with respect to the maximum of the distribution.

5 Applications

Ultimately, the goal of performing frame selection is to try to improve the contrast in the innermost regions, increasing the probability of detecting faint point sources. To quantify the effects of our centering and the frame selection approaches, we studied two targets with known companions that are located at different projected separations from their host star.

thumbnail Fig. 12

S/N of β Pictoris b as a function of the number of principal components used. The solid black line corresponds to the S/N from the data reduction presented in Stolker et al. (2019). The thin solid lines correspond to the full dataset, while the thick solid lines correspond to the same frames used in Stolker et al. (2019) (matched dataset), but both with our different centering options (depending on the panel). Left panel: frames are processed using the star position as the reference center for the full analysis. Middle panel: results when using the AGPM position for the principal component analysis and the star position for the derotation and stacking of the frames. Right panel: center is located at the position of the AGPM for the full analysis.

5.1 Signal-to-noise Ratio Improvements

5.1.1 The Case of β Pictoris b

β Pictoris is a young star that hosts a gas-bearing debris disk (Smith & Terrile 1984; Mouillet et al. 1997; Dent et al. 2014). The star is located at a distance of 19.45 ± 0.05 pc (van Leeuwen 2007) and is a member of the β Pictoris moving group with an estimated age of 22 ± 6Myr (Shkolnik et al. 2017). β Pictoris hosts at least one young giant planet, β Pictoris b, discovered in 2008 via direct imaging (Lagrange et al. 2009, 2010). More recently the presence of another planet, β Pictoris c, has been suggested (Lagrange et al. 2019b). The planet ß Pictoris b is located at 8.900.41+0.23$8.90_{ - 0.41}^{ + 0.23}$ au; it has a period of 20.291.35+0.86$20.29_{ - 1.35}^{ + 0.86}$ yr (Lagrange et al. 2019a) and a mass of 9.32.5+2.6$9.3_{ - 2.5}^{ + 2.6}$ (see Brandt et al. 2021; Snellen & Brown 2018).

Observations with VLΠNaCo and the AGPM were carried out on 2013 February 01 (ESO program ID: 60.A – 9800(J)). They were first presented in Absil et al. (2013), and processed again in Stolker et al. (2019). The summary of the observations is presented in Table 1. It should be noted that the observations were performed under poor weather conditions, which gave us the opportunity to test the effect of frame selection and centering under less than ideal conditions.

The observations were reduced as described in Sect. 2, including our approach for the centering of the science frames. The frame selection is performed following the workflow presented in Sect. 3.3. We tried with different values of ϕ, from +1 to –1 in steps of 0.25. These values, combined with ρ = 0.5, provide us with a wide range of frame selection efficiency, from ~9% to ~65% of the total number of frames (34277). With these cubes, as well as the complete dataset (i.e., without frame selection), we are able to compare our approach with other state-of-the-art post-processing techniques. In this case, we compare our results with those presented in Stolker et al. (2019). We worked with the same frames as in Stolker et al. (2019) to directly compare the effect of the centering on the final images. Stolker et al. (2019) performed a first frame selection where the frames with an unusually high background level were removed from the cube just before aligning the images. Consequently, we used 29 681 frames out of the original 34 277 to compare the effect of different centering techniques. In addition, we also used the entire sequence of 34 277 frames in order to compare the effect of rejecting those frames. In the rest of this work all images are cropped to 45 × 45 pixels.

For each cube of data we measured the S/N of β Pictoris b for a different number of principal components, ranging between 1 and 120 in steps of 1. The signal was measured using an aperture with a radius equal to the FWHM of the PSF measured on the photometric frame, and the noise was estimated using other apertures placed at the same distance to the center of the image. To compute the S/N of the point source, we used the python pipeline module called FalsePositiveModule from the PynPoint package (Amara & Quanz 2012; Stolker et al. 2019), which accounts for low number statistics (Mawet et al. 2014). To better characterize the point-source in every reduced image, we fitted a 2D Gaussian profile to obtain the coordinates of the companion. With these coordinates, weighted by their respective uncertainties, we estimate the mean position to be used in the PynPoint routine.

When performing the principal component analysis we consider three cases for the centering of the science frames. We center them in the position of the AGPM (to later be derotated and stacked) using either the position of the AGPM or of the star as the center of rotation; we also center them at the position of the star for both the principal component analysis and the derotation and stacking processes. The S/N as a function of the number of principal components for the three centering strategies using the full and matched datacubes are presented in Fig. 12. The same figure, but for the frame selection with four different selection criteria, is presented in Fig. 13. It is worth noting that the sky subtraction for the data presented in Stolker et al. (2019) was performed differently than in this work. They estimated the sky contribution using a principal components analysis while we use the sky frames that are closest in time to correct each science frame. Therefore, it is not surprising that the S/N curves show some differences for certain numbers of components used, while the overall trend remains the same.

The S/N estimation does not include an estimate of its uncertainty and the measure of the S/N may vary due to the fluctuation of the noise when removing the stellar PSF and combining the images together. To obtain a first-order estimate of this uncertainty, we first fit a polynomial of degree degree to the values of S/N as a function of the number of components, degree is fixed using the Bayesian information criterion (BIC, Schwarz 1978). The BIC is a criterion for model selection among a finite set of models (in this case a family of polynomials of maximum degree 50), where the increased complexity of the model is penalized (i.e., models with a small number of free parameters are preferred). The model with the lowest BIC value is selected, and then we subtract the resulting polygon from the S/N curve, and measure the standard deviation as an uncertainty estimator from the residuals.

Overall, we find that for β Pictoris b we reach a better S/N (higher values) when the centering used for the principal component analysis is located at the AGPM position, followed by the derotation at the location of the star. This combination of centers for the entire sequence of frames provides the highest S/N, with an improvement of 19.5% ± 4.2% with respect to the analysis presented in Stolker et al. (2019). We note that if we compare the matched dataset (29 681 frames) to this combination of centers, we obtain a slightly higher S/N than presented in Stolker et al. (2019), but within the uncertainties.

Concerning the S/N when applying frame selection, we first note that the S/N agrees with the reduction presented in Stolker et al. (2019) with an increase in the S/N between 40 and 60 principal components. Except for values of ϕ of 1 and 0, and using only the star position in the post-processing analysis, overall we obtain significant improvements in the S/N of β Pictoris b. The best S/N values are 23 ± 1 when using only the star as a reference center, and 22.6 ± 0.6 when using the combination of AGPM and star positions as reference centers. However, considering the uncertainties, we conclude, and in agreement with the full dataset, that the combination of the AGPM and star positions provides the best reduction.

When we set the center for the principal component analysis at the location of the AGPM, we are aligning the torus and the speckles together in all frames. Therefore, the speckles and stellar PSF subtraction become more efficient and result in cleaner individual images. The improvements are most notable in the innermost regions, close to the AGPM, compared to regions that are less affected by speckles and the stellar PSF. This results in higher S/N values, and can also be seen when comparing the final images in Fig. 14. We also note that selecting the location of the AGPM to perform the principal component analysis seems to reduce the self-subtraction effect around the companion, but we do not investigate this further in this study. Our approach to performing the derotation centered at the location of the star is justified by the fact that the rotational center of NaCo, when performing ADI observations, is the central star according to the manual of the instrument6. Derotating at a different location would lead to smearing effects on potential companions when mean- or median-stacking the dataset.

thumbnail Fig. 13

Similar to Fig. 12, but here the thick transparent lines correspond to the reduction with frame selection using ϕ equal to 1 (red), 0.0 (blue), –0.5 (green), and –1.0 (cyan). For all the reductions, ρ = 0.5.

thumbnail Fig. 14

Final images of β Pictoris b using different approaches for the data processing. Left: images of β Pictoris b using the reduction from Stolker et al. (2019), using 49 principal components. Right: same as the left panel, but with our approach for the centering using the entire dataset, using the AGPM position for PCA, and using the star position for the derotation and stacking of the frames. The color scale is linear and is the same for both images.

5.1.2 The Case of RCrAb

R Coronae Australis (R CrA - HIP 93449) is a HAeBe star, one of the brightest in the very young, compact, and obscured Coronet protostar cluster (Taylor & Storey 1984). The star is at an early evolutionary stage and is still embedded in its dusty envelope (Kraus et al. 2009). There are some discrepancies in the estimation of the age, which varies between 0.3 and 3.1 Myr (e.g., Bibo et al. 1992; Forbrich et al. 2006; Meyer & Wilking 2009; Sicilia-Aguilar et al. 2011). The distance to the Corona Australis region is estimated to 154 ± 4 pc in Dzib et al. (2018) using mean trigonometric parallax obtained from Gaia Collaboration (2018). Cugno et al. (2019) and Mesa et al. (2019) identified a close companion in the R CrA envelope: R CrA b. Cugno et al. (2019) reported a mass for the companion in the range 0.10–0.63 M at an angular separation of 196.8 ±4.5 mas (using the first epoch of NaCo L’ observations taken on 2017 May 19), while Mesa et al. (2019) estimated a mass of 0.29 ± 0.08 M with an angular separation of 184 ± 4mas (using SPHERE observations, Kl-K2 dual band imaging, 2018 August 16). In addition, Sissa et al. (2019) reported that the star is a spectroscopic binary, with masses for the primary and secondary of 3 M and 2.3 M, respectively, and a period of ~60 days.

The target was observed as part of the ISPY program, and the observations were presented in Cugno et al. (2019). Here we used the first of the two epochs published in Cugno et al. (2019). It is important to note that the star saturated around the AGPM during the observing sequence (fully saturated at ~ 19 000 counts, see Appendix E). We reduced the observations and estimated the S/N the same way as for β Pictoris b. We used the same central mask to block the central region as in Cugno et al. (2019, 0″05 in radius), and the frames were cropped to the same size (39 × 39 pixels). Even though the target was observed under good weather conditions (see Table 1), we find that for at least half of the time the star is shifted toward a preferential direction with respect to the AGPM, which suggests rather poor centering (see Fig. D.1, and Appendix E for examples of frames with different centering).

We produced cubes similar to those for β Pictoris, and we similarly considered three different centering strategies; using solely the AGPM, solely the star, or both as reference centers (see Section, ring strategy we produced datacubes for the full dataset, or for subsets after performing frame selection. The frame selection was performed using ϕ values between 1 and –1 in steps of 0.25. In the end, a total of 30 different cubes were built, in addition to the cube presented in Cugno et al. (2019) where the authors used a different centering strategy. We produced the final images using a number of principal components between 1 and 40 in steps of 1, for each datacube, since for large numbers of components self-subtraction can become important. Then the S/N of the companion was estimated for all the different reductions in the same way as described for β Pictoirs b in Sect. 5.1.1. In the end, we obtain 31 different S/N curves, one per datacube, and Fig. 15 shows a selection of them for the full dataset and when performing frame selection (keeping ~70, ~55, and ~47% of all frames).

From Fig. 15 we highlight that using the AGPM as the center for the entire process (right panel) yields the worst S/N curve, with a maximum S/N on the order of 10 whether we use frame selection or the full dataset (similar values were obtained by Cugno et al. 2019). In the other two cases (either using the position of the star only, or both the positions of the AGPM and of the star) we note significant improvements on the estimated S/N. For all the cases we find that the S/N is higher in the principal component range of 13–25. In particular, when we use both the AGPM and the star positions in the post-processing analysis, the peak broadens, reaching a maximum S/N of 20.3 when using the full dataset (middle panel of Fig. 15). This broadening of the peak suggests that not only is the noise lower, but also that self-subtraction effects become less significant as the number of principal components increases. In contrast, when using only the star as the center for the entire process, we find that the S/N of the companion is overall lower, suggesting that the PSF subtraction is less adequate in this case despite a narrow peak when using 17 principal components. Given the narrowness of this peak, the result should be taken with caution.

The negative patterns, due to self-subtraction effects, appear more symmetric when we use only the star position or both the star and AGPM positions as the center (see Fig. 16 to compare the images with the best S/N). In Appendix D we show the distribution of relative positions between the star and the AGPM (for all the targets), and as mentioned above the distribution for R CrA is not centro-symmetric with respect to the position of the AGPM. If only the AGPM is used as the center of reference, for instance, the signal of the object might be more diluted in the final image, due to the choice of center for the derotation. For R CrA we conclude that this is the most critical point when post-processing the observations as it can significantly improve the results in the inner regions. Regarding the impact of frame selection, we find that overall the improvements are marginal at best for this specific target.

thumbnail Fig. 15

S/N of R CrA b as a function of the number of components used. The solid black line corresponds to the S/N from the data reduction presented in Cugno et al. (2019). The thin solid lines correspond to the full dataset. The thick transparent lines correspond to the reduction with frame selection using ϕ equal to 1 (orange), 0.5 (pink), and 0.25 (cyan). For all the reductions, ρ = 0.5. In the left panel the frames are processed using the star position as the reference center for the full analysis. The middle panel shows the results when using the AGPM position for the principal component analysis and the star position for the derotation and stacking of the frames. In the right panel the center is located at the position of the AGPM for the full analysis.

5.2 Effects on Asimmetrie Measurements

Astrometric measurements are critical for assessing the co-moving nature of candidate companions. The astrometric position as a function of time allows us to study the orbital motion of giant gaseous planets around their host stars, helping to better constrain their orbital parameters and masses (see, e.g., Wang et al. 2018; Nielsen et al. 2020; Brandt et al. 2021). For this reason, we attempt to quantify here (at least to first order) the effect of the centering and the frame selection on the astrometry of point sources.

The positions of β Pictoris b and R CrA b can be obtained following different approaches. For example, the negative fake companion technique (Lagrange et al. 2010; Marois et al. 2010), which consists of injecting a negative signal at the approximate location of the source of interest in the frame sequence at each parallactic angle. The technique aims to cancel out the companion as much as possible in the final post-processed image and has been used for instance in Stolker et al. (2019) and Cugno et al. (2019) (see Sect. 3.5.2 of Stolker et al. 2019 for further information). There are other methods to determine the astrometric position of companions, for example using the ANgular Differential OptiMal Exoplanet Detection Algorithm (ANDROMEDA, Cantalloube et al. 2015) or the Vortex Imaging Processing (VIP, Gomez Gonzalez et al. 2016), which have been applied in Biller et al. (2021), Jorquera et al. (2021), and Langlois et al. (2021). However, all these methods require considerable computing power, proportional to the size of the field of view and the number of frames. In the case of NaCo we have 10000–30000 frames per target, making these procedures very challenging to carry out. Stolker et al. (2019) and Cugno et al. (2019) addressed this problem by stacking several images together, reducing the number of frames from thousands to hundreds. When using frame selection, we remove frames that are not necessarily consecutive in time, which means we cannot easily stack the images in a similar manner; we therefore opted for the simpler approach of fitting 2D Gaussian profiles to the images.

For each dataset (β Pictoris and R CrA, with or without frame selection, different centering, or original reductions), we computed the final images using a different number of principal components (between 1 and 35 for RCrA, and between 1 and 120 for β Pictoris). For each image we then fitted a 2D Gaussian profile to find the location of the companion (see Appendix F, Figs. F.1 and F.2 for the position of the companion as a function of the number of principal components for β Pictoris b and R CrA b, respectively). The free parameters are the amplitude, center, rotation angle, and width for both axes. For each dataset we then computed the mean value for the X and Y positions, weighted by their corresponding uncertainties. Residual speckles systematically affect the position of the companion in an image using a certain number of components. However, these residuals change while increasing the number of components used, affecting the observed position in a different manner. In this way we can obtain a better approximation of the position of the companion (least affected by this systematic) by calculating the median between all the positions obtained using a wide range of principal components.

Figures 17 and 18 show the results for β Pictoris b and R CrA b where we plot the results for different datasets. Despite small differences, the astrometry remains consistent for all the datasets on each companion (including the dataset from Stolker et al. 2019 and Cugno et al. 2019) at the 3σ level. The figures show that, overall, frame selection reduces in some cases the uncertainties on the astrometric measurements when compared to the reductions for the full datasets. For the astrometry of R CrA b, the dataset from Cugno et al. (2019) agrees better with our solution using only the AGPM. The same is observed when comparing the S/N where both follow the same trend (see Fig. 15). The astrometry uncertainties for R CrA b are larger than for β Pictoris b and, even though all the datasets are consistent with each other, differences of ~ 0.7 pixels are registered corresponding to the magnitude of the asymmetry on the distribution of the differences between star and AGPM positions. Overall, we do not observe significant differences in the final astrometry at the 3σ level. This is the case even without accounting for additional uncertainties, for example on the pixel scale (typically 27.2 ± 0.1 mas pixel−1) or the true north (0.4 ± 0.2 deg). We therefore conclude that the centering approach seems to primarily have an impact for the S/N of the detection, and can reduce the uncertainties associated with the astrometry, but marginally affects the astrometric measurements themselves.

thumbnail Fig. 16

Final images with the highest S/N for R CrA b using different approaches for the data processing. Top row: images using principal components from 13 to 35. Bottom row. same as the top row but for principal components between 1 and 13. Shown are the best S/N for Cugno et al. (2019) (left), for the entire dataset using our centering approach with the combination of the AGPM and star positions (middle), and the frame selection using only the star position as reference center (right). The top and bottom rows have different color scales, but the subpanels in each row are consistent.

6 Summary and Conclusions

In this paper we presented the new CenteR algorithm to improve the performance of coronagraphic NaCo observations at L’ band. One of the fundamental challenges of such observations is to find the location of the star behind the coronagraph (AGPM). Our method makes use of the sky frames to determine the location of the AGPM in the science frames, thus alleviating the degeneracy when modeling the contributions from the central star and AGPM to infer the positions of either. The method described in this paper can be applied to any instrument for which the circular aperture of the coronagraphic mask is visible in all the frames and where the AGPM emits thermally. We also presented an in-depth analysis of the impact of frame selection when dealing with datasets containing tens of thousands of frames. Our approach is based on the modeling results of the star and AGPM and on the pseudo-Zernike moment decomposition of the sky-subtracted frames to select the most homogeneous images. We also investigated the impact of choosing different centers for the principal component analysis and the derota-tion of an angular differential imaging sequence. Our pipeline is tested on two different stars harboring known companions at short and intermediate separations (R CrA and β Pictoris). We also tested our pipeline on a star that hosts a bright circumstellar disk (HD 34282, see Appendix G). However, the analysis of the disk and its detectability and S/N improvements is beyond the objectives of this article.

Using observations of a star with a bright nearby companion, HD 104237, we demonstrated that our approach can successfully recover the true position of the star behind the AGPM, with a ~0.20 pixel dispersion, even if the observing sequence is interrupted. Since our method provides the positions of the AGPM and of the star, we showed that performing the PCA centered at the location of the AGPM followed by the derotation at the location of the star yields better results compared to using only one single position (either the star or the AGPM). Performing PCA at the location of the AGPM removes the diffracted starlight and speckles more efficiently, but the derotation and stacking of the frames should be performed centering at the location of the star to avoid smearing the flux from any nearby point sources. For NaCo the center of rotation coincides with the science target, so that no shifts have to be applied beforehand. We find that for point sources the flux is less smeared, yielding better detections and higher S/N values. Regarding the astrometry, we find that all the centering strategies are consistent with each other within 3σ, while the corresponding uncertainties are smaller when using both the star and the AGPM positions. We compared our reductions with published results, reduced with state-of-the-art pipelines, and found that we are able to improve the S/N of the detection of point sources around ß Pictoris and R CrA by 24 ± 3 and 117 ± 11%. This may also be related to the intrinsic quality of the centering during the observations. Our method yields higher S/N values when the star is poorly centered behind the AGPM and the distribution of the difference between the star and AGPM positions is not centro-symmetric. In contrast, our method reaches similar or slightly better S/N values when the distribution of the difference between the star and the AGPM positions is centro-symmetric.

To the best of our knowledge, this study is the most thorough investigation of the impact of frame selection for datasets consisting of thousands of images. We concluded that frame selection, at best, marginally improves the final reduction. We find that frame selection is most useful when either the observing conditions or the AO correction degraded during the observing sequence. The frames taken under such conditions would dilute the astrophysical signal, and removing them does help decrease the noise, especially in the innermost regions that are more likely to be dominated by speckles. Regarding the astrometry, we find that frame selection can slightly decrease the uncertainties depending on the centering strategy and the percentage of frames kept, but overall provide consistent solutions when compared with the full datasets. On the other hand, in cases for which the centering is poor and anisotropic (i.e., the star moved in a preferential direction during the observing sequence), the centering strategy using the full dataset can have an impact on the astrometry. R CrA is a clear example of this situation, where we estimated a possible bias in the position of about 0.7 pixels or ~ 19 mas when using only the AGPM centering in the entire processing. Nonetheless, for this particular case, the different centering strategies implemented provide results that are consistent with each other within 3σ.

There is a wealth of archival data obtained with an AGPM vector vortex coronagraph. First- and second-generation instruments that use an AGPM, such as VLT/NaCo (Rousset et al. 2003; Lenzen et al. 2003), LBT/LMIRCam (Skrutskie et al. 2010; Defrère et al. 2014) or Keck/NIRC2 (Vargas Catalán et al. 2016; Serabyn et al. 2017) and its previous and ongoing surveys (e.g., LEECH Skemer et al. 2014; Stone et al. 2018), may benefit from a reanalysis using the method described in this study. New generation instruments, VLT/ERIS or ELT/METIS among others, may benefit from the use of an AGPM and the centering method and techniques described in this study. Thanks to the improved centering of the whole dataset, our approach yields better results closer to the AGPM, in the typical range 10–50 au from the star, where the population of giant planets is still poorly constrained.

thumbnail Fig. 17

Position in the final image (45 × 45 pixels) of β Pictoris b. Each subpanel corresponds to a specific dataset (frame selections, full dataset, and matched dataset). The solid red circles correspond to the datasets that use the position of the star as the reference center, while the solid blue squares use only the AGPM, and the solid green triangles use both (star + AGPM). The orange diamond corresponds to the reduction of Stolker et al. (2019). The ellipses correspond to 3σ uncertainties, while the error bars 1σ uncertainties.

thumbnail Fig. 18

Similar to Fig. 17, but for R CrA b (compared with Cugno et al. 2019).

Acknowledgements

The authors thank the anonymous referee for a thorough review. The comments helped improving the manuscript, especially with respect to the astrometric measurements, but also the overall structure of the paper. N.G., J.O. and A.B. acknowledge financial support by ANID, – Millennium Science Initiative Program – NCN19_171. N.G. acknowledges grant support from project CONICYT-PFCHA/Doctorado Nacional/2017 folio 21170650. J.O. acknowledges financial support from the Universidad de Valparaíso, and from FONDE-CYT Regular (grant 1180395). A.B. acknowledges support from FONDECYT Regular 1190748. A.M. and A.Q. acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (MU 4172/1-1). Th.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. G.M.K is supported by the Royal Society as a Royal Society University Research Fellow. T.S. acknowledges the support from the ETH Zurich Postdoctoral Fellowship Program. S.P.Q. and G.C. thank the Swiss National Science Foundation for financial support undergrant number 200021_169131.

Appendix A Motion of the Circular Aperture as a Function of Time

Figure A.1 shows the motion of the center of the circular aperture as a function of time for all four targets. The science frames (colored solid circles) and the sky frames (colored solid triangles) move on the same random path as a function of time. Even with that random path the motion of the circular aperture in consecutive science-sky frames is less than 0.1 pixels.

thumbnail Fig. A.1

Motion in time of the circular aperture for the four stars studied in this work. Left panel: Position of the center of the circular aperture along the X-axis as a function of time for all four targets. If the sequence was interrupted (because of an instrumental problem or zenith avoidance region), each sequence is plotted in a different panel. The error bars correspond to the 1σ uncertainties. The colored solid circles correspond to science frames and the colored solid triangles to the sky frames. Right panel: Same as the left panel, but for the Y-axis.

Appendix B AGPM: Circular Aperture Motion as a Function of Time

Figure B.1 shows the difference between the AGPM position and the circular aperture center as a function of time for the X- and Y-axes. The relation shows small deviations of higher polynomial order than a constant behavior, but which are within the data uncertainties at the 2σ level. To avoid overfltting, we chose the simplest solution and adopted a constant value calculated from the mean weighted by the respective uncertainties. In addition, Figure B.1 shows the calculated value for every target that provides the direct transformation between the center of the circular aperture and the AGPM position. In the case of HD 104237 we calculated two different transformations for each of the observing sequences (before and after the interruption). We further note that as the center of the AGPM is located near the center of the circular aperture, a possible rotation of the circular aperture in time (instead of a linear shift) will not have a strong impact on the differential motion between the circular aperture and the AGPM.

thumbnail Fig. B.1

Differences between the position of the AGPM and the center of the circular aperture as a function of time (sky frames only), for the X- and Y-axis (left and right, respectively), for all four targets. The error bars shown represent the 1er uncertainties. The solid horizontal lines correspond to the mean weighted by the uncertainties, and the colored areas correspond to the 1er uncertainty.

Appendix C The Effect of H in Frame Selection

The homogeneity parameter H is obtained using the reconstruction of pseudo-Zernike moments for the nonradial and azimuthal components (see Eq. 6 and Fig. 6). The reconstruction is a linear combination where the base corresponds to the pseudo-Zernike moments. For this reason, it is possible to consider the square root sum of the coefficients as a constant (e.g., normalizing the root square sum to 1). When the star is well centered behind the AGPM, the torus formed is mostly homogeneous, where the radial component dominates (larger values for |m| = 0 coefficients and therefore for TorusS/N, see Eq. 7 and Fig. 6 left panel). When the star is misaligned with respect to the AGPM, the azimuthal components dominate in the reconstruction (higher values for |m| = 1 and 2 coefficients, lower H values, see Eq. 6 and Fig. 6 right panel). Then, H and TorusS/N are correlated with each other and both provide similar information. Figure C.1 shows the normalized TorusS/N as a function of the star position with respect to the AGPM center along the X-axis in the left panel, and in the right panel the values for H with the same horizontal axis. The red dots correspond to the selected frames (~ 55%), the gray dots to the rejected ones (~ 45%), and the blue dots are the frames that, having been selected in the original frame selection (red dots), are excluded when we also consider the H parameter for the selection, using a cutoff value below 0.6. We note that both selection criteria are mutually consistent, with a 3.6% discrepancy (corresponding to ~ 2% of all the frames). The location of these blue dots in both graphics shows that both TorusS/N and H classify those frames as low-quality frames with respect to the other ones. Therefore, the incorporation of the H parameter does not have a strong impact on data processing since it provides information similar to and consistent with TorusS/N.

thumbnail Fig. C.1

Example of registration and selection of frames with and without considering the H parameter for HD 34282. Left: TorusS/N as a function of the difference between the star position and AGPM center along the X axis. Right: H asa function of the difference between the star position and AGPM center along the X axis. The red dots correspond to the frame selection presented in this work, using ϕ = 0.5 and ρ = 0.5 corresponding to ~ 55% of the frames, while the gray dots are the rejected frames. The blue dots are the additional rejected frames when H is also considered in the frame selection criterion. The number of blue dots corresponds to 3.6% and 2% of the first selection (red points) and the full dataset, respectively.

Appendix D Relative Positions and Torus S/N

The left panel of Figure D.1 shows the position of the stars with respect to the AGPM for all the sources. The right panel shows the de-normalized TorusS/N for a more direct comparison between all the observations at different weather conditions. In particular, R CrA shows a distribution that is not centro-symmetric. This implies that the choice of the center becomes more critical, affecting the S/N and the astrometric measurements of the companion. In contrast, β Pictoris shows a more uniform distribution (azimuthally homogeneous) where we do not observe significant differences in the results for the S/N and positions of the companion.

thumbnail Fig. D.1

Distribution of the frame registration and density distribution of the relative positions for each star. Left panel: Distribution of the difference between the positions of the star and of the AGPM for all four stars included in this study. The ΔX and ΔY correspond to ~ 1.5 times the median absolute deviation. Right panel: TorusS/N for all the targets in a common reference frame. The contours contain 10%, 30%, 50%, 68%, and 95% of the data.

Appendix E Frames of R CrA

Regions of the frames of R CrA are saturated, especially close to the AGPM in the entire observing sequence, resulting in a stronger signal of the torus and speckles around the AGPM. Figure E.1 shows the distribution of the star position with respect to the center of the AGPM along both axes (top left panel) and with respect to the TorusS/N and homogeneity H (top right and bottom left panel, respectively). Examples of the different images at different separations are shown in the bottom right panel, corresponding to examples from good to poor centering (red triangles, green squares, and cyan circles, respectively). In this particular case, both TorusS/N and H are less informative in distinguishing between frames with good and poor centering. However, the distribution of the star position with respect to the AGPM is a more robust criterion to select the most homogeneous and well-centered frames.

thumbnail Fig. E.1

Frame registration for R CrA for six of the parameters used in the frame selection, and examples of frames with different centering. Top left: Distribution of the star position with respect to the AGPM center along both axes. Top right: TorusS/N as a function of the star position with respect to the AGPM center along the X-axis. Bottom left: Homogeneity H as a function of the star position with respect to the AGPM center along the X-axis Bottom right: Examples of images with good (red triangles), intermediate (green squares), and poor (cyan circles) centering. The yellow stars correspond to the star position, while the red crossed circles give the AGPM position. The images have the same color scale.

Appendix F Position of Companions as a Function of Number of Components

We made a 2D Gaussian fitting to each of the companions (β Pictoris b and R CrA b) to know their approximate positions in each of the post-processed images according to the number of components used. Since this method of knowing the position of the object in question is very sensitive to residual speckles and then to the number of components used, we took all the images to calculate the weighted average by the uncertainties associated with each position. Figure F.1 shows an example for β Pictoris b of the results of our fitting for the X- and Y-axis (solid points with their error bars), as well as the calculated average position and its associated uncertainty (horizontal solid line and the colored area, respectively), for the case of full datacubes. Figure F.2 shows an example for R CrA b in the same format. These figures show the trend of the positions of each of the objects as a function of the number of components. It is noticeable that as the number of components increases, the associated uncertainty also increases due to the fact that there is less flux available for the fitting (effect of self-subtraction).

thumbnail Fig. F.1

Each panel shows the position for X (top panels) and Y (bottom panels) using our 2D Gaussian fitting for β Pictoris b. From left to right, the center used for our reduced data correspond to: only the star, the star and the AGPM, only the AGPM, and at the end the dataset from Stolker et al. (2019). The solid horizontal line corresponds to the weighted mean by the respective uncertainties (the vertical line in each solid point), with the respective associated uncertainty (colored horizontal area).

thumbnail Fig. F.2

Same as Fig. F.1, but for R CrA b.

Appendix G Extended sources

The advantage of using direct imaging, compared with Polarimetrie observations, is that it is possible to obtain the full intensity of extended structures. However, it is usually a challenge for angular differential imaging as it removes part of the flux and can introduce dark regions at the edges. In most cases the circumstellar disks are significantly removed in the final image especially for disk with inclinations lower than 60° (see Milli et al. 2012) as a disk with low inclination will be incorporated in the reference PSF to be subtracted. Different efforts have been put into addressing and solving this problem, either from the observations themselves (e.g., RDI, Smith & Terrile 1984) or in post-processing techniques (see, e.g., LOCI, Lafrenière et al. 2007). More recently, Ren et al. (2018, 2020) implement a new algorithm, the data imputation using sequential non-negative matrix factorization (DI-sNMF), which is a model-free solution in a post-detection disk scenario. Broadly speaking, the idea is to transform the signal (disk or point sources) into a missing data problem, assigning to that region the PSF signal (following the statistical properties). Then it should be possible to remove both the stellar PSF and the speckles in every image, avoiding problems of self-subtraction. However, none of these studies investigate the effect of the centering or frame selection on the properties of the resulting disk image (signal, structure, self-subtraction), and detectability (e.g., noise distribution, significance of the detection).

We study here the impact of both the centering and frame selection when the target shows extended emission in its vicinity. From the previous sections, we know that the centering has a significant impact on the results for point sources, but the effect of frame selection is marginal. However, frame selection may improve the speckle suppression in the case of disks. If we keep the most homogeneous, well-centered frames, we should be able to decrease the number of components needed, overall decreasing the contribution of residual speckles in the inner regions. Furthermore, by removing the frames with a poor quality (poor centering or weather conditions) the contrast should improve in the inner regions where the disk is located. We tested our reduction strategy on the HAeBe star HD 34282. The star is located at a distance of 359 ± 5pc (Gaia Collaboration 2018), with an age of 6.412.58+1.92$6.41_{ - 2.58}^{ + 1.92}$ Myr (Merín et al. 2004). The disk was first resolved with ALMA observations (van der Plas et al. 2017), and first imaged at infrared wavelengths by Launhardt et al. (2020) as part of the ISPY program using VLΠNaCo, and by de Boer et al. (2021) as part of the DISK program in Polarimetrie mode at near-infrared using VLT/SPHERE. The reduction is performed in the same way as for the previous targets in this paper. The details related to the observations are summarized in Table 1. We performed the frame selection using ϕ and ρ as 0.5, keeping ~ 50 % of the 32 825 available frames.

Figure G.1 shows reduction examples for three different centering configurations (top to bottom), using 18 principal components. From left to right, we show the final images when using the full dataset, with our frame selection, and a control reduction, selecting at random the same number of frames as for the middle column. Quantifying the S/N for extended emission is more difficult than for a point source; therefore, we just compare the images visually, paying special attention to the noise level and to self-subtraction effects. For instance, when using the full dataset the dark pattern caused by self-subtraction is quite significant, along with spots toward the east. In contrast, the reduction using frame selection shows a better removal of several speckles, consequently the structure of the disk is overall better recovered. This can be explained by the presence of a significant number of poorly centered frames in the entire dataset, which increases the noise in the innermost regions and dilutes the faint signal from the disk.

When performing frame selection we may remove a significant number of frames, and we therefore introduce gaps in the sequence of frames. These gaps correlate with the observing condition or the AO performance. As a consequence, we most likely remove series of consecutive frames, which may act as a protection angle7 (see LOCI, Lafrenière et al. 2007), therefore decreasing the effect of self-subtraction. To understand whether the result presented in Figured (middle panel) is due to the frame selection or the result of a nonhomogeneous sampling in time, we carried out the following test. First, from the cube with frame selection (using ϕ = 0.5 and ρ = 0.5), we counted the number and sizes of the gaps introduced in the sequence. We then cloned this distribution to obtain a new sample of gaps, which we could randomly introduce in the entire dataset. This new subsample had the same number of frames as the previous one. The result is shown in Figured (right column). Overall, the signal from the disk appears fainter, the western side seems to be more disconnected from the trace of the disk, and the speckles on the eastern side are still visible, while they had disappeared when using our frame selection approach. This result strongly suggests that our strategy to perform frame selection does improve the final reduction, and that it is not solely related to the sampling of frames as a function of time. We also note that, in agreement with the previous section, using only the AGPM position as the reference center does not yield the best final image. Without a proper quantification of the S/N of the disk, it remains difficult to estimate which of the other centering strategies provides the best reduction.

thumbnail Fig. G.1

Comparison of the final images for the disk around HD 34282 using 18 principal components for the full dataset, frame selection, and random selection of frames. Left panel: Final images using the entire dataset, but implementing three different centers: the AGPM position only, both the AGPM and the star positions, and only the star position. Middle panel: Same as the previous panel, but with frame selection, using ~ 50% of the data. Right panel: Same as the middle panel but with a random selection of frames (see text for details). Red arrows indicate the most noticeable lingering shadows and speckles, as well as the most notable disk features. The yellow contours correspond to the continuous image of the observations made with ALMA (van der Plas et al. 2017), and correspond to 10 and 20 mJy levels.

References

  1. Absil, O., Milli, J., Mawet, D., et al. 2013, A&A, 559, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Absil, O., Mawet, D., Delacroix, C., et al. 2014, SPIE Conf. Ser., 9148, 91480M [NASA ADS] [Google Scholar]
  3. Absil, O., Mawet, D., Karlsson, M., et al. 2016, SPIE Conf. Ser., 9908, 99080Q [NASA ADS] [Google Scholar]
  4. Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [Google Scholar]
  5. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bhatia, A. B., & Wolf, E. 1954, Math. Proc. Cambridge Philos. Soc., 50, 40 [CrossRef] [Google Scholar]
  7. Bibo, E. A., The, P. S., & Dawanas, D. N. 1992, A&A, 260, 293 [NASA ADS] [Google Scholar]
  8. Biller, B. A., Apai, D., Bonnefoy, M., et al. 2021, MNRAS, 503, 743 [NASA ADS] [CrossRef] [Google Scholar]
  9. Böhm, T., Catala, C., Balona, L., & Carter, B. 2004, A&A, 427, 907 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2019, A&A, 624, A87 [EDP Sciences] [Google Scholar]
  11. Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, ApJ, 898, L16 [Google Scholar]
  12. Bowler, B. P. 2016, PASP, 128, 102001 [Google Scholar]
  13. Brandt, G. M., Brandt, T. D., Dupuy, T. J., Li, Y., & Michalik, D. 2021, AJ, 161, 179 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, A&A, 582, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Chauvin, G. 2018, SPIE Conf. Ser., 10703, 1070305 [NASA ADS] [Google Scholar]
  16. Cugno, G., Quanz, S. P., Launhardt, R., et al. 2019, A&A, 624, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. de Boer, J., Ginski, C., Chauvin, G., et al. 2021, A&A, 649, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Defrère, D., Absil, O., Hinz, P., et al. 2014, SPIE Conf. Ser., 9148, 91483X [Google Scholar]
  19. Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014, Science, 343, 1490 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dzib, S. A., Loinard, L., Ortiz-León, G. N., Rodríguez, L. F., & Galli, P. A. B. 2018, ApJ, 867, 151 [Google Scholar]
  21. Feigelson, E. D., Lawson, W. A., & Garmire, G. P. 2003, ApJ, 599, 1207 [Google Scholar]
  22. Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020, A&A, 637, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Forbrich, J., Preibisch, T., & Menten, K. M. 2006, A&A, 446, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Fusco, T., Rousset, G., Sauvage, J. F., et al. 2006, Opt. Express, 14, 7515 [NASA ADS] [CrossRef] [Google Scholar]
  26. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [Google Scholar]
  27. Garcia, P. J. V., Benisty, M., Dougados, C., et al. 2013, MNRAS, 430, 1839 [NASA ADS] [CrossRef] [Google Scholar]
  28. Gerard, B. L., & Marois, C. 2016, SPIE Conf. Ser., 9909, 990958 [Google Scholar]
  29. Gomez Gonzalez, C. A., Absil, O., Absil, P. A., et al. 2016, A&A, 589, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Grady, C. A., Woodgate, B., Torres, C. A. O., et al. 2004, ApJ, 608, 809 [NASA ADS] [CrossRef] [Google Scholar]
  31. Groff, T. D., Kasdin, N. J., Limbach, M. A., et al. 2015, SPIE Conf. Ser., 9605, 96051C [NASA ADS] [Google Scholar]
  32. Hagelberg, J., Ségransan, D., Udry, S., & Wildi, F. 2016, MNRAS, 455, 2178 [NASA ADS] [CrossRef] [Google Scholar]
  33. Huby, E., Baudoz, P., Mawet, D., & Absil, O. 2015, A&A, 584, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Huby, E., Bottom, M., Femenia, B., et al. 2017, A&A, 600, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Hunziker, S., Quanz, S. P., Amara, A., & Meyer, M. R. 2018, A&A, 611, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Jorquera, S., Pérez, L. M., Chauvin, G., et al. 2021, AJ, 161, 146 [NASA ADS] [CrossRef] [Google Scholar]
  37. Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890 [NASA ADS] [CrossRef] [Google Scholar]
  38. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Kraus, S., Hofmann, K. H., Malbet, F., et al. 2009, A&A, 508, 787 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770 [Google Scholar]
  41. Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148 [CrossRef] [Google Scholar]
  42. Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [Google Scholar]
  44. Lagrange, A. M., Boccaletti, A., Langlois, M., et al. 2019a, A&A, 621, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019b, Nat. Astron., 3, 1135 [NASA ADS] [CrossRef] [Google Scholar]
  46. Langlois, M., Gratton, R., Lagrange, A. M., et al. 2021, A&A, 651, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Launhardt, R., Henning, T., Quirrenbach, A., et al. 2020, A&A, 635, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Lenzen, R., Hartung, M., Brandner, W., et al. 2003, SPIE Conf. Ser., 4841, 944 [Google Scholar]
  49. Ma, H., Qiao, Y., & Shen, C. 2017, SPIE Conf. Ser., 10463, 1046316 [NASA ADS] [Google Scholar]
  50. Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci., 111, 12661 [NASA ADS] [CrossRef] [Google Scholar]
  51. Males, J. R., Close, L. M., Miller, K., et al. 2018, SPIE Conf. Ser., 10703, 1070309 [NASA ADS] [Google Scholar]
  52. Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541 [Google Scholar]
  53. Marois, C., Doyon, R., Racine, R., & Nadeau, D. 2000, PASP, 112, 91 [Google Scholar]
  54. Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
  55. Marois, C., Macintosh, B., & Véran, J.-P. 2010, SPIE Conf. Ser., 7736, 77361J [Google Scholar]
  56. Mawet, D., Absil, O., Delacroix, C., et al. 2013, A&A, 552, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [Google Scholar]
  58. Merín, B., Montesinos, B., Eiroa, C., et al. 2004, A&A, 419, 301 [Google Scholar]
  59. Mesa, D., Bonnefoy, M., Gratton, R., et al. 2019, A&A, 624, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Meyer, M. R., & Wilking, B. A. 2009, PASP, 121, 350 [NASA ADS] [CrossRef] [Google Scholar]
  61. Milli, J., Mouillet, D., Lagrange, A. M., et al. 2012, A&A, 545, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896 [Google Scholar]
  63. Mugnier, L. M., Cornia, A., Sauvage, J.-F., et al. 2009, J. Opt. Soc. Am. A, 26, 1326 [NASA ADS] [CrossRef] [Google Scholar]
  64. Murphy, S. J., Lawson, W. A., & Bessell, M. S. 2013, MNRAS, 435, 1325 [Google Scholar]
  65. Nielsen, E. L., De Rosa, R. J., Wang, J. J., et al. 2020, AJ, 159, 71 [Google Scholar]
  66. Noll, R. J. 1976, J. Opt. Soc. Am., 66, 207 [Google Scholar]
  67. Pueyo, L. 2016, ApJ, 824, 117 [Google Scholar]
  68. R Core Team 2019, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria [Google Scholar]
  69. Rajwa, B., Dundar, M., Irvine, A., & Dang, T. 2013, IM: Orthogonal Moment Analysis, r package version 1.0 [Google Scholar]
  70. Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104 [Google Scholar]
  71. Ren, B., Pueyo, L., Chen, C., et al. 2020, ApJ, 892, 74 [Google Scholar]
  72. Rousset, G., Lacombe, F., Puget, P., et al. 2003, SPIE Conf. Ser., 4839, 140 [Google Scholar]
  73. Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118 [NASA ADS] [CrossRef] [Google Scholar]
  74. Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Samland, M., Bouwman, J., Hogg, D. W., et al. 2021, A&A, 646, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Sauvage, J.-F., Fusco, T., Petit, C., et al. 2016, J. Astron. Teles. Instrum. Syste., 2, 025003 [NASA ADS] [CrossRef] [Google Scholar]
  77. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  78. Serabyn, E., Huby, E., Matthews, K., et al. 2017, AJ, 153, 43 [NASA ADS] [CrossRef] [Google Scholar]
  79. Shkolnik, E. L., Allers, K. N., Kraus, A. L., Liu, M. C., & Flagg, L. 2017, AJ, 154, 69 [Google Scholar]
  80. Sicilia-Aguilar, A., Henning, T., Kainulainen, J., & Roccatagliata, V. 2011, ApJ, 736, 137 [NASA ADS] [CrossRef] [Google Scholar]
  81. Sissa, E., Gratton, R., Alcalà, J. M., et al. 2019, A&A, 630, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Skemer, A. J., Hinz, P., Esposito, S., et al. 2014, SPIE Conf. Ser., 9148, 91480L [NASA ADS] [Google Scholar]
  83. Skrutskie, M. F., Jones, T., Hinz, P., et al. 2010, SPIE Conf. Ser., 7735, 77353H [NASA ADS] [Google Scholar]
  84. Smith, B. A., & Terrile, R. J. 1984, Science, 226, 1421 [Google Scholar]
  85. Snellen, I. A. G., & Brown, A. G. A. 2018, Nat. Astron., 2, 883 [Google Scholar]
  86. Soummer, R., Hagan, J. B., Pueyo, L., et al. 2011, ApJ, 741, 55 [Google Scholar]
  87. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [Google Scholar]
  88. Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543 [Google Scholar]
  89. Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Stone, J. M., Skemer, A. J., Hinz, P. M., et al. 2018, AJ, 156, 286 [Google Scholar]
  91. Taylor, K. N. R., & Storey, J. W. V. 1984, MNRAS, 209, 5P [CrossRef] [Google Scholar]
  92. Teh, & Chin, R. T. 1988, IEEE Trans. Pattern Anal. Mach. Intell., 10, 496 [CrossRef] [Google Scholar]
  93. van der Plas, G., Ménard, F., Canovas, H., et al. 2017, A&A, 607, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. van Leeuwen, F. 2007, A&A, 474, 653 [CrossRef] [EDP Sciences] [Google Scholar]
  95. Vargas Catalán, E., Huby, E., Forsberg, P., et al. 2016, A&A, 595, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72 [EDP Sciences] [Google Scholar]
  97. Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192 [Google Scholar]
  98. Xuan, W. J., Mawet, D., Ngo, H., et al. 2018, AJ, 156, 156 [NASA ADS] [CrossRef] [Google Scholar]
  99. Zernike, F. 1934, MNRAS, 94, 377 [Google Scholar]

1

According to Appendix A.12, page 116 in the VLT/SPHERE User Manual Issue: 106.

2

PynPoint version 0.6.2: https://pynpoint.readthedocs.io/en/latest/

4

Or 10 times σimage below the median when the observations are performed with windowing. In this case, there are many fewer pixels outside the aperture, and therefore the median is no longer representative of the background level.

5

We consider the previous fit done using Eq. (1) and their uncertainties to constrain the σagpm.

6

According to Sect. 5.8, page 64 in the VLT/NACO User Manual Issue: 102.

7

The protection angle consists of removing some of the adjacent frames in the time series. For the ith frame in the sequence, only the frames for which the paralactic rotation is larger than ∆θ are used to build the reference, where ∆θ depends on the angular size of the PSF at a specific separation from the center. The selected frames are used in the ADI, and this approach minimizes the effect of self-subtraction.

All Tables

Table 1

Log of the observations.

Table 2

Astrometric measurements of HD 104237 B.

All Figures

thumbnail Fig. 1

Example of a sky image showing the circular aperture in blue (observing date 2016 November 08. The top right inset shows a zoom-in on the central part, with the AGPM thermal emission in the center. The white elongated spots correspond to dust or dirt in the optical path of the instrument or on the camera, and are emitting at near-IR wavelengths.

In the text
thumbnail Fig. 2

Examples of what the AGPM and the circular aperture look like in a sky image. Left: modeling of the thermal emission profile of the AGPM along the X- and Y-axes to obtain the coordinate center for a cleaned sky image. Right: modeling of the circular aperture along the X- and Y-axes to estimate the coordinates of its center in the same image. The colored area corresponds to a 95% confidence interval.

In the text
thumbnail Fig. 3

Position of the center of the circular aperture (with respect to the camera) during the observation sequence for sky and science images (red and blue, respectively), for HD 34282 (left), β Pictoris (middle), and R CrA (right); the axis scales are 1:5, 1:1, and 1:3, respectively. The grid represent pixels, and the colored crosses correspond to the uncertainty obtained in the fit.

In the text
thumbnail Fig. 4

Position of the center of the circular aperture (blue triangles) and AGPM (red circles) with respect to their median value, for sky frames only. Left: position for HD 34282. Middle: same, but for β Pictoris. Right: same, but for R CrA. The axis scales of the left, middle, and right panels are 1:8, 1:1, and 1:3, respectively. The colored crosses correspond to the uncertainty obtained on the fit.

In the text
thumbnail Fig. 5

Modeling example of sky-subtracted science images for HD 34282 and β Pictoris using both the negative and positive 2D Gaussian. For each sub-panel, the observations are shown on the left, while the best-fit model is shown on the middle, and the residuals in the right. The homogeneity (H) and the S/N of the torus (TorusS/N) is labeled for each frame. The color scale is the same for each subpanel. The white cross for each center panel marks the location of the AGPM.

In the text
thumbnail Fig. 6

Two examples of reconstructed images using the pseudo-Zernike moments. For each panel (made up of six subpanels): top left: observations; top center: reconstructed image using m up to 12; top right: residuals of the modeling; bottom left: reconstruction of the image using only m = 0, related to radial and uniform azimuthal distribution; bottom center: reconstruction of the image using |m| = 1 and 2, corresponding to the adopted inhomogeneous contribution; bottom right: reconstruction of the image using 1 ≤ |m| ≤ 12, related to all the azimuthal contributions. For each subpanel the color scale is the same.

In the text
thumbnail Fig. 7

Examples of the distribution and registration of frames as a function of different parameters. Left: S/N of the torus as a function of the difference between the positions of the star and of the AGPM along the X-axis. Right: relative position of the star with respect to the position of the AGPM. For both panels the star is HD 34282 and the color-coding corresponds to the goodness of fit. The crosses in both panels correspond to the typical uncertainties. The contours contain 93 and 68% of the data in the left panel, and 93, 68, and 50% of the data in the right panel.

In the text
thumbnail Fig. 8

Examples of selected or rejected frames after applying our selection criteria. Left: normalized TorusS/N as a function of the difference between the positions of the star and the AGPM along the X-axis. The color-coding corresponds to the frame selected using ϕ = 0.5 and ρ = 0.5 (red dots) and rejected frames (gray dots). The orange circles correspond to three examples of frames selected and the cyan triangles to examples of rejected frames. The small panel at the top of shows the histogram of both distributions in logarithmic scale. Right: three examples of rejected frames (top row) and selected frames (bottom row). The red crossed circles correspond to the position of the AGPM, while the yellow star gives the derived position of the star obtained from the negative and positive 2D Gaussian fitting. The color scale is the same for all six images.

In the text
thumbnail Fig. 9

Percentage of selected frames as a function of the free parameter ϕ for HD 34282. The different lines show the behavior for three different values of ρ.

In the text
thumbnail Fig. 10

Density contours of the position of the close companion after derotating each frame by the parallactic angle and the position angle of the companion, and after subtracting its angular separation. The derotation is either done at the position of the star inferred from our modeling (left) or at the position of the AGPM (right). The black circle corresponds to the uncertainties on the angular separation and the position angle of the companion at the peak of the distribution. The grid shown in each panel represent pixels, and the colored crosses correspond to the typical uncertainty. The contour lines in each panel enclose the areas containing 99.7, 95.5, 68.5, 50, 30, and 10% of the sample with respect to the maximum of the distribution.

In the text
thumbnail Fig. 11

Comparison of the star positions derived from the close companion and those obtained by our fitting approach (negative and positive 2D Gaussian fitting labeled NP2DG in the figure). The small top and side panels show the kernel density estimators (KDE) for the X- and Y-axes, respectively. The grid represent pixels, and the colored cross corresponds to the typical uncertainty. The contour lines enclose the areas containing 95.5, 68.5, 50, 30, 10, and 1% of the sample with respect to the maximum of the distribution.

In the text
thumbnail Fig. 12

S/N of β Pictoris b as a function of the number of principal components used. The solid black line corresponds to the S/N from the data reduction presented in Stolker et al. (2019). The thin solid lines correspond to the full dataset, while the thick solid lines correspond to the same frames used in Stolker et al. (2019) (matched dataset), but both with our different centering options (depending on the panel). Left panel: frames are processed using the star position as the reference center for the full analysis. Middle panel: results when using the AGPM position for the principal component analysis and the star position for the derotation and stacking of the frames. Right panel: center is located at the position of the AGPM for the full analysis.

In the text
thumbnail Fig. 13

Similar to Fig. 12, but here the thick transparent lines correspond to the reduction with frame selection using ϕ equal to 1 (red), 0.0 (blue), –0.5 (green), and –1.0 (cyan). For all the reductions, ρ = 0.5.

In the text
thumbnail Fig. 14

Final images of β Pictoris b using different approaches for the data processing. Left: images of β Pictoris b using the reduction from Stolker et al. (2019), using 49 principal components. Right: same as the left panel, but with our approach for the centering using the entire dataset, using the AGPM position for PCA, and using the star position for the derotation and stacking of the frames. The color scale is linear and is the same for both images.

In the text
thumbnail Fig. 15

S/N of R CrA b as a function of the number of components used. The solid black line corresponds to the S/N from the data reduction presented in Cugno et al. (2019). The thin solid lines correspond to the full dataset. The thick transparent lines correspond to the reduction with frame selection using ϕ equal to 1 (orange), 0.5 (pink), and 0.25 (cyan). For all the reductions, ρ = 0.5. In the left panel the frames are processed using the star position as the reference center for the full analysis. The middle panel shows the results when using the AGPM position for the principal component analysis and the star position for the derotation and stacking of the frames. In the right panel the center is located at the position of the AGPM for the full analysis.

In the text
thumbnail Fig. 16

Final images with the highest S/N for R CrA b using different approaches for the data processing. Top row: images using principal components from 13 to 35. Bottom row. same as the top row but for principal components between 1 and 13. Shown are the best S/N for Cugno et al. (2019) (left), for the entire dataset using our centering approach with the combination of the AGPM and star positions (middle), and the frame selection using only the star position as reference center (right). The top and bottom rows have different color scales, but the subpanels in each row are consistent.

In the text
thumbnail Fig. 17

Position in the final image (45 × 45 pixels) of β Pictoris b. Each subpanel corresponds to a specific dataset (frame selections, full dataset, and matched dataset). The solid red circles correspond to the datasets that use the position of the star as the reference center, while the solid blue squares use only the AGPM, and the solid green triangles use both (star + AGPM). The orange diamond corresponds to the reduction of Stolker et al. (2019). The ellipses correspond to 3σ uncertainties, while the error bars 1σ uncertainties.

In the text
thumbnail Fig. 18

Similar to Fig. 17, but for R CrA b (compared with Cugno et al. 2019).

In the text
thumbnail Fig. A.1

Motion in time of the circular aperture for the four stars studied in this work. Left panel: Position of the center of the circular aperture along the X-axis as a function of time for all four targets. If the sequence was interrupted (because of an instrumental problem or zenith avoidance region), each sequence is plotted in a different panel. The error bars correspond to the 1σ uncertainties. The colored solid circles correspond to science frames and the colored solid triangles to the sky frames. Right panel: Same as the left panel, but for the Y-axis.

In the text
thumbnail Fig. B.1

Differences between the position of the AGPM and the center of the circular aperture as a function of time (sky frames only), for the X- and Y-axis (left and right, respectively), for all four targets. The error bars shown represent the 1er uncertainties. The solid horizontal lines correspond to the mean weighted by the uncertainties, and the colored areas correspond to the 1er uncertainty.

In the text
thumbnail Fig. C.1

Example of registration and selection of frames with and without considering the H parameter for HD 34282. Left: TorusS/N as a function of the difference between the star position and AGPM center along the X axis. Right: H asa function of the difference between the star position and AGPM center along the X axis. The red dots correspond to the frame selection presented in this work, using ϕ = 0.5 and ρ = 0.5 corresponding to ~ 55% of the frames, while the gray dots are the rejected frames. The blue dots are the additional rejected frames when H is also considered in the frame selection criterion. The number of blue dots corresponds to 3.6% and 2% of the first selection (red points) and the full dataset, respectively.

In the text
thumbnail Fig. D.1

Distribution of the frame registration and density distribution of the relative positions for each star. Left panel: Distribution of the difference between the positions of the star and of the AGPM for all four stars included in this study. The ΔX and ΔY correspond to ~ 1.5 times the median absolute deviation. Right panel: TorusS/N for all the targets in a common reference frame. The contours contain 10%, 30%, 50%, 68%, and 95% of the data.

In the text
thumbnail Fig. E.1

Frame registration for R CrA for six of the parameters used in the frame selection, and examples of frames with different centering. Top left: Distribution of the star position with respect to the AGPM center along both axes. Top right: TorusS/N as a function of the star position with respect to the AGPM center along the X-axis. Bottom left: Homogeneity H as a function of the star position with respect to the AGPM center along the X-axis Bottom right: Examples of images with good (red triangles), intermediate (green squares), and poor (cyan circles) centering. The yellow stars correspond to the star position, while the red crossed circles give the AGPM position. The images have the same color scale.

In the text
thumbnail Fig. F.1

Each panel shows the position for X (top panels) and Y (bottom panels) using our 2D Gaussian fitting for β Pictoris b. From left to right, the center used for our reduced data correspond to: only the star, the star and the AGPM, only the AGPM, and at the end the dataset from Stolker et al. (2019). The solid horizontal line corresponds to the weighted mean by the respective uncertainties (the vertical line in each solid point), with the respective associated uncertainty (colored horizontal area).

In the text
thumbnail Fig. F.2

Same as Fig. F.1, but for R CrA b.

In the text
thumbnail Fig. G.1

Comparison of the final images for the disk around HD 34282 using 18 principal components for the full dataset, frame selection, and random selection of frames. Left panel: Final images using the entire dataset, but implementing three different centers: the AGPM position only, both the AGPM and the star positions, and only the star position. Middle panel: Same as the previous panel, but with frame selection, using ~ 50% of the data. Right panel: Same as the middle panel but with a random selection of frames (see text for details). Red arrows indicate the most noticeable lingering shadows and speckles, as well as the most notable disk features. The yellow contours correspond to the continuous image of the observations made with ALMA (van der Plas et al. 2017), and correspond to 10 and 20 mJy levels.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.