© H. Le Coroller et al. 2020

1 Introduction

Most of the 4100 exoplanets detected to date have been found using indirect methods, such as the radial velocity technique and photometric transits. It is indeed extremely difficult to detect planet light that is drowned in the much brighter diffracted light from its host star. Jupiter and Earth like planets are about 10⁸ to 10¹⁰ fainter than their parent star in the visible band.

However, owing to a combination of eXtreme Adaptive Optics (ExAO), innovative coronagraphs, differential imaging, and sophisticated post-processing algorithms, direct imaging instruments have been able to detect and characterize young giant planets at large separation (>~10 au). Across two decades of exoplanetary science in direct imaging, dozens of dedicated surveys have been carried out around young, nearby stars (Chauvin 2018b). These have led to the discovery of the first planetary mass companions in the early 2000s at large distances of ≥100 au. The implementation of differential techniques, starting in 2005, enabled the breakthrough discoveries of closer and lighter planetary mass companions such as HR 8799 bcde (10, 10, 10, and 7 M_Jup at 14, 24, 38, and 68 au, respectively; Marois et al. 2008, 2010), β Pictoris b (8 M_Jup at 8 au; Lagrange et al. 2009), Fomalhaut b (< 1 M_Jup at 177 au; Kalas et al. 2008; still debated), HD 95086 b (5 M_Jup at 52 au; Rameau et al. 2013).

The current generation of ExAO planet imagers (Macintosh et al. 2014; Jovanovic et al. 2015; Beuzit et al. 2019), Gemini Planet Imager (GPI), Subaru Coronagraphic Extreme Adaptive Optics (SCExAO), Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE) are now equipped with integral field spectrographs offering exquisite near-infrared spectra of young giant planets to unveil the physical processes at play in their atmospheres and a link to their mechanisms of formation. For the first time, these instruments have reached a contrast level of ≈ 10⁻⁵ at a separation of about 200 to 900 mas, enabling the detection of the new young planets 51 Eri b (2 M_Jup at 13 au; Macintosh et al. 2015), HIP 65426 b (9 M_Jup at 92 au; Chauvin et al. 2017b), and PDS 70 b (9 M_Jup at 29 au; Keppler et al. 2018).

An important finding from these high-contrast imaging surveys in recent years has been the low occurrence rate of giant planets beyond 30 au (${0.6}_{-0.5}^{+0.7}\%$; see Bowler 2016). Today, the Gemini Planet Imager Exoplanet Survey (GPIES) and the SPHERE High-contrast ImagiNg survey for Exoplanets (SHINE) large surveys of about 600 observed stars indicate that this scarcity extends down to 10 au (Nielsen et al. 2019; Vigan et al. 2020), suggesting that the bulk of the giant planet population is typically located between 1 au and 10 au.

A prime goal of the future surveys will also be to bridge the gap with indirect techniques by imaging young Jupiters down to the snowline at about 3–5 au, depending on stellar type. The next generation of instruments such as SPHERE+ (Boccaletti et al. 2020) aim to reach contrasts of at least 10⁻⁵ at ≈ 100 mas, which represents an improvement of a factor of 3 in terms of angular separation with respect to SPHERE (Beuzit et al. 2019). K-Stacker, together with the SPHERE+ capability, has the potential to achieve the core of the Jupiter-mass planets population at ≈ 3 au (i.e., ≈ 10⁻⁶ at ≈ 60 mas) in 10− 100 h of exposure time, by providing anadditional factor of 3−10 in contrast. In terms of observing strategy, for planets at less than 10 au from a star at 10 pc to 20 pc, the orbital motion becomes comparable to the width of a 10 m telescope diffraction-limited Point Spread Function (PSF) in about 30 days. Taking into account observing constraints and weather statistics because higher contrast can only be reached during the best nights with a seeing <0.6″, multiple observations of very interesting young, nearby systems are usually spread over several days, months, and years (see case of β Pictoris; Lagrange et al. 2019a). In this case, the orbital motion of the potential planets makes a simple co-addition of the different images suboptimal in terms of pure detection, if not impossible.

The Keplerian motion of the planet has to be taken into account during the combination. On the Extremely Large Telescopes (ELTs), the situation will be even worse. The Point Spread Function (PSF) will be four to five times smaller than with the 10 m telescope class such as Very Large Telescope (VLT), Keck, Gemini, Subaru, Large Binocular Telescope (LBT), etc. These ELTs will be used to search for planets at very small separations (below 10 au) with first-light instruments (Chauvin 2018a). In this case, the Keplerian motion could become non-negligeable in a matter of only a few days (Males et al. 2013).

Following similar principles previously applied to the search of new moons in solar systems such as Hippocamp, the seventh innermost moon of Neptune (Showalter et al. 2019), the Keplerian-Stacker (K-Stacker) algorithm described by Le Coroller et al. (2015) is an observing strategy and method of data reduction applied to nearby stars. This technique consists in combining high-contrast images recorded during different nights, accounting for the orbital motion of the putative planet that we are looking for. Even if an individual image does not reveal the planet, we showed that an optimization algorithm like K-Stacker can be used to properly align the images according to Keplerian motion; for instance, 10–50 images taken over the course of several months or years. The resulting gain in signal-to-noise ratio (S/N) can allow for the detection of planets that are otherwise unreachable. This method can be used in addition to angular differential imaging (ADI) and spectral differential imaging (SDI) techniques (Racine et al. 1999; Marois et al. 2006) or any other high-contrast data reduction method designed to further improve the global detection limit. As a byproduct of the optimization algorithm, K-Stacker also directly provides the orbital parameters of the detected planets.

Ultimately, the main goal of K-Stacker would be to directly drive the observing strategy and scheduling of current and future planet imagers in which exposures would be split over several nights to maximize the detection performances. In Paper I (Nowak et al. 2018), using simulated VLT-SPHERE observations, we have shown that when the total number n of available images is large enough to get $\sqrt{n}\times$ (S/N) ≥ 7 (where (S/N) is the signal-to-noise levels in individual frames), the K-Stacker algorithm is able to detect the planet with a high level of reliability >90%. The number of false positives were low but the simulated images did not reproduce instrumental speckle noise and angular spectral differential imaging (ASDI) reductions. The main goal of this paper is to validate the K-Stacker algorithm on real data obtained on sky reduced by the most recent ADI and ASDI Principal Component Analysis (PCA) algorithms.

In Sect. 2, we describe the observations used in this paper, which come from the SPHERE SHINE survey (Chauvin et al. 2017a). In Sect. 3, we present the results of a blind test in which K-Stacker was used to search for fake planets hidden in real InfraRed Dual-band Imager and Spectrograph (IRDIS) observations reduced by a PCA ADI algorithm. We also study the capability of K-Stacker to recover the correct orbital parameters despite the typical errors encountered in real data, such as instrumental true north offsets or stellar mass and distance uncertainties. In Sect. 4, we show that K-Stacker is able to recover the known companions β Pictoris b and HD 95086 b and show that the orbital parameters retrieved by K-Stacker are in agreement with the values found in the literature. We present in Sect. 5 the first K-Stacker searches for new planets around HD 95086 and β Pictoris, two targets which have been repeatedly observed during the SHINE survey. In Sect. 6, we discuss the strategy of future K-Stacker observations. Section 7 gives our final conclusions.

2 Description of the observations used in this paper

All the observations used in this paper come from the SHINE survey done with the Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE) instrument at the focal plane of the VLT-UT3 (Beuzit et al. 2019). The SHINE program (Chauvin et al. 2017b) is a very high-contrast near-infrared survey of more than 600 young, nearby starsaimed at searching for and characterizing new planetary systems. The goal of this project is also to find statistical constraints on the rate, mass, and orbital distributions of the giant planet population at large orbits. Evenif the SHINE observations were not organized for the K-Stacker “philosophy”, in which we plan to observe fewer stars but over more epochs, the number of observed targets reduced homogeneously by the SPHERE Data Center (Delorme et al. 2017; Galicher et al. 2018) allow us to test K-Stacker for the first time in real conditions.

The SPHERE instrument includes an extreme adaptive optics system, several types of coronagraphs, and three subsystems: IRDIS, the Integral Field Spectrograph (IFS), and the Zurich Imaging POLarimeter (ZIMPOL). In this paper we use observationscoming from IFS and IRDIS, which were designed to cover the near-infrared range for an efficient search of young planets (Beuzit et al. 2019). In the SHINE survey, two Apodized Lyot Coronagraphs (ALC) configurations were used: APO1-ALC2 (N-ALC-YJH-S) coronagraph with a focal plane mask of a diameter of 185 mas and the APO1-ALC3 (N-ALC-Ks) coronagraph with a focal plane mask of a diameter of 240 mas optimized for the IRDIFS and IRDIFS-EXT modes, respectively (Beuzit et al. 2019). In these modes, IRDIS and IFS work simultaneously (Claudi et al. 2008; Zurlo et al. 2014) and IRDIS is used in a dual band imaging configuration (Dohlen et al. 2008; Vigan et al. 2010).

All the IRDIS observations used in the blind test described in Sect. 3 were acquired using the APO1-ALC2 coronagraph. In Sects. 4 and 5, we also use the observations on two emblematic stars HD95086 and β Pictoris that have been observed regularly in the SHINE survey using the APO1-ALC2 and ALC3 coronagraph (see Table B.1) to constrain the orbital parameters of the known planets b around these targets.

3 K-Stacker computation on real SHINE data

3.1 Setup of the blind test

In order to extend the demonstration of the K-Stacker algorithm on simulated IRDIS datasets presented in Paper I, we performed a new analysis on real IRDIS observations obtained during the SHINE survey. The methodology followed during this new blind experiment is similar to the approach discussed in Paper I, with the exception of an additional PCA ADI reduction step. To stay as close as possible to the conditions of Paper I and to demonstrate the true potential of K-Stacker, we injected planets in fake “K-Stacker runs”, each made of ten observations taken from the SHINE survey. As most stars have not been observed that many times during the survey, we created the K-Stacker runs by combining observations acquired on different but similar targets. In particular, we were careful in selecting observations of stars with similar magnitudes, and acquired in similar conditions such as seeing and adaptive optics (AO) performances. We created a total of five such fake K-Stacker runs, using a total of 50 different images from the survey.

Fake planets were then injected in the raw observations of each run by doing the following:

• 1.

  Randomly draw one run in which no planet is injected.

• 2.

  For each of the four other runs, draw a set of random orbital parameters (see Table 1 for an overview of orbital laws used), and inject the planet in each raw observation of the run according to the orbit drawn. The target star is assumed to have a mass of 1 M_⊙, and to be located at 10 pc. The planet is injected at a random contrast uniformly drawn in the range [5 × 10⁻⁶, 4 × 10⁻⁷].

• 3.

  Reduce each observation of each run using the PCA ADI tools of the SPHERE Data Center (Delorme et al. 2017; Galicher et al. 2018).

The process was repeated six times to create a total of 30 fake K-Stacker runs, in which 6 have no planets, and 24 have a planet injected randomly. The resulting runs were then classified based on the perfectly recombined (S/N)_tot of the injected fake planets. We note that the S/N definition when using K-Stacker can be confusing, as we need to distinguish three different values: the S/N in each individual ADI observation of the run (simply S/N), the S/N after the K-Stacker recombination (S/N)_KS, and the optimal S/N that could be achieved if the orbit was perfectly known, and the images perfectly recombined (S/N)_tot. The 30 runs were then divided into the following three groups:

• Group I: composed of 13 runs for which ${(S/N)}_{\text{tot}}<2$ for the injected fake planets (or non injected fake planets). These planets can be considered as undetectable, that is, equivalent to no planet injected.

• Group II: composed of 9 runs for which $2<{(S/N)}_{\text{tot}}<12$ for the injected fake planets (low-S/N regime)

• Group III: composed of 8 runs for which ${(S/N)}_{\text{tot}}>12$ for the injected fake planets and therefore easily detectable.

3.2 Blind check

Following the principles of Paper I, after running the K-Stacker algorithm, we asked an independent observer, who was not aware of the presence or not of fake planets in the images, the injection phase, the S/N, or the orbital parameters, to check the final K-stackerrecombined solutions and to assign for each run one of the three flags: “no detection”, “planet candidate”, and “possible candidate, more observations required”.

Among the 13 sets of Group I (i.e., ${(S/N)}_{\text{tot}}<2$), 8 were correctly flagged by the observer as no detection (i.e., true negatives); the five other solutions were flagged as “possible candidates, more observations needed”. Because of the nature of the blind test performed, in which the same limited number of high-contrast observations were used multiple times, all the runs are not fully independent, and these five cases correspond to two truly independent candidates emerging from the noise. The typical K-Stacker S/N for these five cases was (S∕N)_KS ≈ 6, and although the observer did not claim a detection from these runs, in reality, they would probably have led to more observation time being spent on the targets. We consider these two cases as false positives.

Among the nine runs of Group II (i.e., $2<{(S/N)}_{\text{tot}}<12$), five solutions were flagged as true detections by the observer and were indeed true positives. Two solutions at ${(S/N)}_{\text{tot}}=3.3$ and ${(S/N)}_{\text{tot}}=3.6$ were found at (S∕N)_KS ≈ 6, and flagged by the observer as possible detections, with more observations required. These two solutions correspond to the same false positives as described above. Two other solutions at ${(S/N)}_{\text{tot}}=3$ and ${(S/N)}_{\text{tot}}=4$ were flagged as no detection by the observer, although the K-Stacker (S/N)_KS was found to be between 7 and 8. In these two last cases, the orbits found by K-Stacker pass well at less than one pixel of the true positions of the fake planet in three images (epochs), but miss the planet along the other epochs. The contribution from the fake planet to the (S/N)_KS explains the relatively high values obtained. We count these ambiguous cases as false negatives. Finally, all eight runs of Group III (i.e., ${(S/N)}_{\text{tot}}>12$) were flaggedby the observer as detections and are indeed true positives.

In Fig. 1, we report the histogram of the planets found and missed from Groups I and II. For clarity, solutions of Groups III with ${(S/N)}_{\text{tot}}>12$ are not shown. All these solutions with ${(S/N)}_{\text{tot}}>12$ correspond to true positives. Even if the statistics are poorer than in Paper I, we still find that K-Stacker is able to recover all the planets with ${(S/N)}_{\text{tot}}>6.5$ i.e., S∕N ≃ 2 in individual observations. With only five truly independent runs of 10 observations used to create the 30 fake K-Stacker runs, it is difficult to draw solid conclusions regarding the false positive rate. One of the two false alarms of Fig. 1 was found to be very close to the coronagraphic mask, where the false positive probability may be significantly higher (see also Sect. 5). However, in every case, the observer did not claim a detection from the available observations, but merely suggested that the K-Stacker solution was possibly a planet and requested more observations.

Fig. 1 Distribution of the planet candidates found and missed as a function of the (S/N)tot given a perfect recombination of the images. The true negatives, all grouped at S∕N = 0, are not shown.

3.3 Extracted orbital parameters

The result of the K-Stacker algorithm is a list of orbits sorted by S/N (see Paper I for detail). The S/N of the ≈ 50−100 first orbits is a relatively flat function before decreasing by steps (see example in Fig. 2). For ourstudy of the orbital parameters, we only kept the orbits before this drop in (S/N)_KS (i.e., 68 first orbits in theexample of Fig. 2). We checked that for the true positives of the SHINE blind test the first K-Stacker solutions before the decreaseof (S/N)_KS always contained a set of orbital parameters that pass well by the orbit of injection (see Fig. 4). We show the orbital parameters in 2D maps (Fig. 4) to be able to compare the K-Stacker solutions with the Markov chain Monte Carlo (MCMC) technique on the positions (Beust et al. 2016).

We also give a “mean” solution with its standard deviation (black cross in Fig. 4). For all the planets detected by K-Stacker in the SHINE blind test (true positives in green of Fig. 1), the mean solution of the orbital parameters is at maximum three standard deviations from the parameters of injection. This mean orbit always passes at less than ≈ 0.6 pixels from the true positions of the fake planet in the images (see Table A.1). Figure 3 shows also that the K-Stacker (S/N)_KS found with the mean solution of the orbital parameters is equal within the error bars to the (S/N)_tot of a perfect recombination. To conclude, K-Stacker has recovered the orbital parameters well, although the planet has traveled over a maximum of15−38% of its total orbital period (see Fig. 4 and Table A.1).

Fig. 2 K-Stacker (S/N)KS derivative infunction of the orbit number sorted by S/N.

Fig. 3 K-Stacker (S/N)KS divided by (S/N)tot of a perfect recombination as a function of (S/N)tot.

3.4 K-Stacker tolerances on the errors of real data

To limit the computation time, K-Stacker does not include true north error, or stellar mass and distance as free parameters. These values are fixed values set by the user. In this section, we investigate the capability of K-Stacker to converge despite possible errors on these parameters.

3.4.1 Tolerance on the stellar mass error and impact on the orbital parameters

To test the robustness of K-Stacker to errors on the stellar mass, we performed an experiment, in which we kept the mass of the target star to a fixed value M = 1 M_⊙ when calculating the orbit of the injected planets, but changed the mass used by K-Stacker to M + δM. For a circular face-on orbit of semimajor axis a, the orbital velocity is given by $$V_{\text{orb}}=\sqrt{\frac{GM}{a}}.$$(1)

Thus, a small error δM on the mass M of the star directly translates to an error δV_orb on the velocity given by $$\delta V_{\text{orb}}=\frac{1}{2}\sqrt{\frac{G}{a}}\frac{\delta M}{M^{1/2}}.$$(2)

This error on the orbital velocity leads to a K-Stacker estimate of the position in each image that “drifts” with time, compared to the real position of the planet. Assuming that the algorithm can always play on appropriate parameters to minimize the total error by properly aligning the center of the time series, we can calculate the typical mean position error ΔX_mean for a sequence of N observations acquired at a regular time interval ΔT over a total period T = (N − 1)ΔT as follows: $$\Delta X_{\text{mean}}=\frac{1}{N}{\sum }_{k=1}^{k=N}k\Delta T\delta V_{\text{orb}}=\frac{T}{2}\delta V_{\text{orb}}.$$(3)

For a star at a distance d, the associated angular error is then $$\Delta {\theta }_{\text{mean}}=\frac{1}{4}\frac{T}{d}\sqrt{\frac{G}{a}}\frac{\delta M}{M^{1/2}},$$(4)

where d = 10 pc, T = 3 yr, a ≃ 5 au, and M = 1 M_⊙, an error of 0.1 M_⊙ on the mass of the star leads to a mean angular error of 21 mas, similar to the size of the SPHERE PSF. For most of the stars of the SHINE survey the mass is known with an accuracy better than 10% (Desidera et al. 2015), so the effect should remain limited. Nevertheless, to study the impact of this error on the performance of K-Stacker, we used the algorithm on a fake run while explicitly adding an error on the stellar mass. In Fig. 5, we give the mean distance between the solution found by K-Stacker and the injected position of the planet, as a function of the error on the stellar mass δM. Interestingly, K-Stacker is able to tolerate a relatively large error on the stellar mass (− 0.2 < δ M < 0.2), which should in principle lead to an error of more than a PSF on the position of the planet (see Eq. (4)). This could be explained by the fact that the algorithm somehow compensates for the wrong mass by altering some of the orbital parameters. Figure 6, which gives the retrieved orbital parameters as a function of δM, indicates that both a and e are strongly correlated with δM, and thus that varying these parameters can indeed help to compensate for the error on the stellar mass.

Although this has not been studied in depth in this work, Fig. 5 indicates that K-Stacker can potentially be used to retrieve the mass of the central star together with the other parameters. The maximum S/N value of (S∕N)_KS = 6.34 is obtained at δM = 0 and, in the region of small mass errors (δM∕M < 10%), the (S/N)_KS function shows a clear drop of up to 0.2 to 0.4 when departing from δM = 0, with no other apparent local maximum. The situation is more complicated for higher initial error on the mass, with a distinct secondary maximum of (S/N)_KS around δM = 0.2 M_⊙ visible in Fig. 5. The reason for the existence of such a secondary maximum is not clear and the grid sampling could play an important role in this situation. Furthermore, it remains to be determined whether this apparent local maximum exists in the full parameter space, or if it is an effect of the projection of the (S/N)_KS function against δM only. A gradient descent taking into account all the orbital parameters and the stellar mass simultaneously could potentially avoid this apparent secondary maximum. But this remains to be demonstrated, and a significant reworking of the algorithm is necessary to constrain the stellar mass on targets for which it is largely uncertain. This is out of scope of the current paper but could be considered, if necessary, for example in the case of low-mass stars for which the mass is not always well known.

Fig. 4 Example of orbital parameters coming from one true positive solution of the SHINE blind test (see Sect. 3.2). In each subplot, each point corresponds to the parameters of one orbit found by K-Stacker with the color of its (S/N)KS value. In the diagonal, the histograms for each orbital parameter are shown. The dark cross corresponds to the mean of the K-Stacker orbits. The size of the dark cross is 1σ of each orbital parameter. The star is at the position of the orbital parameters used to inject the fake planet.

Fig. 5 Results of the K-Stacker algorithm when including an error on the stellar mass used by the algorithm to calculate the orbital positions at each epoch. The dots (resp. dashes) show the mean (resp. maximum)distance between the true position and the position found by K-Stacker in all the images as a function of the error on the mass. These dots are color coded to indicate the corresponding (S/N)KS. The black dashed line indicates the size of the FWHM of the instrumental PSF (i.e., K-Stacker has missed the planet at least in one image when a small line is above the dashed line).

3.4.2 Tolerance on the true north error

True north (TN) gives the absolute rotational orientation of observations. An error on the TN is responsible for a rotational misalignment between the different images used by K-Stacker, resulting in an apparent deviation of the planet motion from Keplerian orbits. For reference, the typical true north error in SPHERE, estimated using a well-defined observing strategy of a given astrometric field (Maire et al. 2016), does not exceed 0.1 deg. For NaCo, similar values of 0.1−0.2 deg were obtained over more than a decade (Chauvin et al. 2012, 2015). To check the consequences of a TN error on theK-Stacker algorithm, we simulated several observations taking into account the real distribution of the TN errors measured on SPHERE (see Table A.2). The maximum TN error is 0.08 deg, with a standard deviation of 0.026 deg, which corresponds to 0.4 mas, or ≈1∕100 of the full width at half maximum (FWHM) of the PSF at the edge of the area corrected for AO. At that level, this error should not affect the performance of K-Stacker in any noticeable way. Indeed, we find that K-Stacker has no problem recovering the planets, and gives the same final S/N and orbital parameters, even when multiplying the SPHERE TN errors by a factor of 10.

3.4.3 Tolerance on the stellar distance error

The stellar distances are known from the HIPPARCOS (van Leeuwen 2007) or Gaia missions (Luri et al. 2018). The typical error on the parallax given by Gaia Data Release 2 (Gaia DR2) for bright sources (mag < 14) is δϖ < 0.1 mas (Luri et al. 2018). For a star at 10 pc (ϖ = 100 mas), this translates to a relative error on the distance of δd∕d = δϖ∕ϖ = 0.1%.

To determine the position of the possible planet in each image, K-Stacker first calculates the position of the planet around its central star and then projects this position on the detector. The projected position is directly proportional to d⁻¹. Given that the corrected field for AO on SPHERE is typically 1″, the projection error induced by a 0.1% error on the distance of the star can be at most 1 mas (i.e., much smaller than the instrumental PSF). Consequently, this error is negligible for K-Stacker.

Fig. 6 Difference between the orbital parameters found by K-Stacker and the real values of injection, as a function of the error on the stellar mass.

Fig. 7 Best recombined image resulting from the K-Stacker run on β Pictoris. At each epoch, the images are rotated and shifted to put the planet on its periastron position found by K-Stacker, and the frames are co-added. The planet b is detected at a (S/N)KS level of 24.5.

4 Applying K-Stacker to known exoplanets

In this section, we present the first results obtained with K-Stacker on two real planets: β Pic b and HD 95086 b. Each of these emblematic objects has been observed multiple times with SPHERE during the SHINE survey, and together they provide a good test of K-Stacker in different conditions. β Pic b is on an edge-on orbit with a significant orbital motion. HD 95086 b is rather on pole-on configuration and has moved by only a few PSFs over the five years of IRDIS monitoring.

4.1 β Pic b

Eleven IFS observations of β Pictoris, spread over more than three years between 2015 and 2018, were available in the SHINE survey (Table B.1). We reduced these data with a PCA ASDI algorithm (Mesa et al. 2015). Although in this case, with a mean S/N of 7.42, the planet is clearly detected in each individual image, K-Stacker was set up to look blindly for planets in the range of orbital parameters given in Table 2.

The planet β Pic b is detected at a total K-Stacker recombined (S∕N)_KS = 24.5 (see Fig. 7), a gain of a factor 3.3 compared to the individual PCA ASDI reduced observations. For a set of 11 observations, this gain is optimal.

Figure 8 gives the distribution of the 89 best orbits found by K-Stacker in the parameter space. The black crosses on the different subplots give the position of the mean value of these 89 orbits, with the associated 1σ spread (see also Table 3). For comparison, the red cross gives the best estimates and 1σ uncertainties from Lagrange et al. (2019a), converted to the reference system used in K-Stacker.

Since K-Stacker does not implement a proper MCMC exploration of the parameter space, the statistical meaning of the corner plots presented in Fig. 8 is not straightforward. But these pseudo-corner plots share some interesting similarities with the results of a more classical approach to fitting, as presented in Lagrange et al. (2019a): a clear V-shaped correlation between the semi-major axis a and the eccentricity e, related to a degeneracy on the position of periastron/apoastron; a well constrained edge-on inclination; and an eccentricity distribution that peaks at e < 0.1. Overall, the orbital solution resulting from the K-Stacker run is in good agreement with the recent results of Lagrange et al. (2019a). The only significant difference is on t₀ found by Lagrange et al. (2019a), which is near the apoastron of the K-Stacker mean solution (see red and dark crosses for t₀ in Fig. 8). This ambiguity in the periastron/apoastron is reinforced by the small eccentricity and the incomplete coverage of the orbit. This will be solved with further monitoring.

4.2 HD 95086 b

For HD 95086, a total of eight observations are available in the SHINE survey between May 2015 and May 2019, and all of these were obtained in good conditions in the H part of the spectrum with IFS (see Table B.1). The data were reduced using a PCA ASDI algorithm.

The algorithm was set to search for planets in the range of parameters given in Table 2. The planet HD 95086 b is detected (Fig. 9) at a recombined (S∕N)_KS = 9.97. The mean S/N in the individual PCA ASDI reduced observations was 3.67. Thus, the S/N is improved by a factor 2.72 by K-Stacker. For a series of eight observations, this is again very close to the optimal case. The mean orbital solution found by K-Stacker is presented in Table 3, and the associated pseudo-corner plot can be found in Fig. 10. Although 1.4% of the orbital period of HD 95086 b has been covered by the observations, the orbital parameters are relatively well constrained. Again, the pseudo-corner plots of K-Stacker share some similarities with the results of a more classical approach to fitting, as presented in Chauvin et al. (2018). Within the error bars, the mean values of the orbital parameters found by K-Stacker (see dark and red crosses of Fig. 10) are equal to the solutions coming from the MCMC technique described in Chauvin et al. (2018).

Fig. 8 Histograms and 2D diagrams of the β Pictoris b orbital parameters found by K-Stacker. Left, top, and bottom: scale of the orbital parameters. Right: scale of the histograms. 89 points in each 2D diagram corresponding to the 89 orbits with the highest (S/N)KS found by K-Stacker. The color of each point gives the (S/N)KS indicated at right. The dark cross indicates the mean value of the orbital parameters with their error bars. The red cross shows the higher probability density found by an MCMC technique in Lagrange et al. (2019a) converted in the K-Stacker referential. The origin of the t0 K-Stacker date is the 01∕12∕2014.

Fig. 9 Best recombined image resulting from the K-Stacker run on HD 95086. At each epoch, the images are rotated and shifted to put the planet on its periastron position found by K-Stacker, and the frames are co-added. The planet b is detected at a (S/N)KS level of 9.97.

5 Searching for inner exoplanets

In this Section, we focus on the search for additional inner planets around HD 95086 and β Pictoris using the SHINE reduced data (Table B.1). Indeed, the presence of one or two additional inner giant planets is suspected considering the double-belt architecture of HD 95086 (Chauvin et al. 2018) and that a ~ 9 M_Jup planet β Pic c has been found using radial velocity, at 2.7 au from the star (Lagrange et al. 2019b). To search for additional companions in these systems, we proceed using the same K-Stacker algorithm, in which we introduce an initial step of masking theknown imaged planet b from the individual images. We search in the inner area of HD 95086 and β Pictoris with the parameters given in the Table 4.

For both HD 95086 and β Pictoris, we find a bright secondary feature, close to the coronagraphic mask (see Fig. 11). For HD 95086, the feature is located at a = 16.7 au and has a low (S∕N)_KS = 4.4. Its spread shape indicates a probable false positive, but a more detailled study will be done in a future work (Desgrange et al., in prep.) to try to determine if it is a true detection or not.

In the case of β Pictoris, the compact feature is found at (S∕N)_KS = 4.8 (Fig. 11) with the orbital parameters given at Table 5. This solution is compatible with the orbital parameters inferred from the following radial velocities (Lagrange et al. 2019b): a = 2.6948 au, e = 0.243, t₀ = −7.87206 yr, θ₀ = 4.62 rad (where t₀ and θ₀ were expressed in the K-Stacker referential).

The semimajor axis found by K-Stacker is larger by 0.3 au but the PSF of the planet c is partially masked at several epochs by the coronagraphic mask and can therefore disturb the solution of K-Stacker. However, the parameters Ω =−2.0 rad ± 1.5 rad and i = +1.85 rad ± 0.86 rad give an orbit that is misaligned with the disk that is very difficult to explain from a dynamical point of view. In the blind test (Sect. 3.2), we had two ambiguous cases in which the orbits found by K-Stacker passed well by the planet positions only at three to four epochs over 10. Thus, even if the orbit found by K-Stacker is not correct, a part of the light of this “bright” structure (Fig. 11) could come from β Pictoris c. In all these cases, the probability of a true detection at (S∕N)_KS ≤ 5 is smaller than 50%.

Fig. 10 Histograms and 2D diagrams of the HD 95086 b orbital parameters found by K-Stacker. Left, top, and bottom: scale of the orbital parameters. Right: scale of the histograms. The 129 points in each 2D diagram correspond to the 129 orbits of the higher (S/N)KS found by K-Stacker. The color of each point gives the (S/N)KS indicated at right. The dark cross indicates the mean value of the orbital parameters with their error bars. The red cross shows the higher probability density found by an MCMC technique in Chauvin et al. (2018) converted in the K-Stacker referential. The origin of the t0 K-Stacker date is 05∕05∕2015.

6 Strategy of future K-Stacker observations

A fundamental question for any potential future surveys with K-Stacker will be to select the optimum number of observations per target. The blind tests realized on simulated IRDIS images (Nowak et al. 2018) and on real data (see Sect. 3) show that an (S∕N)_KS > 7 must be reached to claim a true detection with high confidence (i.e., small rate of false alarm). An S∕N > 7 is also required to be able to do a precise spectral analysis and confirm a true detection by characterizing the physicalproperties of the planet (e.g., atmosphere composition, temperature, surface gravity, and clouds detections); see for example the characterization of the planets around the stars HR 8799, β Pic, HD 95086, 51Eri, HIP 65426, PDS 70 (Bonnefoy et al. 2016; Chilcote et al. 2017; De Rosa et al. 2016; Chauvin et al. 2017b; Müller et al. 2018). But, further work is required to develop a K-Stacker exoplanet statistical analysis tool: using a similar approach as with the multipurpose exoplanet simulation system (MESS) algorithm (Bonavita et al. 2012), we should be able to inject numerous fake planets in series of observations to compute a probability of detection in function of the planet mass by using model mass-luminosity relationships (Baraffe et al. 2003). This approach will allow us to give an upper mass limit for planets hidden in the data even when K-Stacker does not detect anything.

The total exposure time required for K-Stacker also depends on the number of observations already done. For example, to confirm the detection of β Pictoris c that has been found at (S∕N)_KS = 4.8 in 11 exposures, a minimum of $11*[{(7/4.8)}^2-1]\approx 12$ new observations are needed to reach (S∕N)_KS ≈ 7. ; this number of exposures could be slightly reduced by using radial velocity constraints to schedule the β Pictoris c observations when the planet is fully out of the coronagraphic mask.

Figure 3 and Table A.1 shows that K-Stacker is able to recombine in a close to optimal way a series of observations with very different orbital parameters. Thus, in principle, the total S/N of a K-Stacker run should only depend on the total exposure time, and not on the number of individual exposures in which it is divided. In these conditions, better constraints on the orbital parameters could be obtained at fixed total exposure time by taking as many exposures over a period as long as possible. But, to reach higher contrast (i.e., per individual exposure and in the final K-Stacker recombined image), a minimum exposure time by observation is required to have enough field-of-view rotation for a good ADI substraction; that is, in each observation, a parallactic rotation of at least one PSF FWHM at the position of the planet must be reached to subtract the instrumental speckles with ADI efficiently. The trade-off comes from the difficulty in setting up the observation strategy, telescope overheads, and star declination (to adjust the minimum exposure of individual observations).

For future imaging surveys with K-Stacker on SPHERE+ (Boccaletti et al. 2020) and/or the ELTs instruments, a software to optimize observing schedules will have to be developed. Inspired from the algorithm used in the SHINE survey (Lagrange et al. 2016), the goal of this K-Stacker scheduler will be to compute the minimum exposure time required at each epoch for each star to have an efficient ADI subtraction, but at the same time to split the observations over at least 10− 20% of the orbital period of the searched planets to get accurate orbital parameters.

Fig. 11 Recombined images obtained with K-Stacker, when searching for additional companions around β Pictoris and HD 95086, two systems observed in the SHINE survey. At each epoch, the images are rotated and shifted to put the detection on its periastron position found by K-Stacker, and the frames are co-added. In each case, a bright spot can be seen near the coronagraphic mask (red arrow), which corresponds to the best solution found by the algorithm. In the case of β Pic, the corresponding (S/N)KS is 4.8. For HD 95086, the (S/N)KS is 4.4.

7 Conclusions

For the first time, we tested the K-Stacker algorithm on real IRDIS and IFS observations, reduced with the PCA ADI and ASDI algorithms, in a similar fashion as in the SHINE survey. From a blind experiment, in which planets are injected on random orbits in the images before the PCA ADI reduction and recovered by the K-Stacker algorithm, we conclude that the detection statistics are similar to what was previously obtained with simulated non-ADI reduced images: the success rate is close to 100% when searching for companions for which the recombined (S∕N)_KS ≃ 9, and it drops significantly at (S∕N)_KS ≃ 5.

Using data on two targets repeatedly observed during the SHINE survey (β Pic, observed 11 times; HD 95086, observed 8 times), we have shown that the K-Stacker algorithm is good at recovering the known companions β Pic b and HD 95086 b. K-Stacker also provides orbital parameter estimates in agreement with current literature values. The gain provided by K-Stacker recombination is very close to the square root of the number of observations combined (i.e., close to optimal).

We also searched for additional substellar companions around these two stars. We found two bright features, corresponding to two possibleplanets orbiting within the orbit of the known companions. However, these two features are found at low (S/N)_KS where it is difficult to reach a conclusion concerning a true detection or not. A dedicated analysis will be done on the peak found around HD 95086 in a future work (Desgrange et al., in prep.). The c candidate detected by K-Stacker around β Pictoris without using prior information from the radial velocity detection of β Pictoris c is on a trajectory compatible with the orbital parameters found by Lagrange et al. (2019b). But, the K-Stacker orbit is misaligned with the disk and the (S∕N)_KS ≈ 5 is not high enough to claim that it is a true detection. Despite the relatively large error bars on the euler angles (Ω, i, θ₀) the 1− σ agreement between the putative K-stacker detection and the radial velocity solution is encouraging and more observations are required to constrain the orbital parameters and to determine if at least a part of the light of this detection comes from β Pictoris c.

Acknowledgements

SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF–Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). We also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). It has also been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). M.R.M., H.M.S., and S.D. are pleased to acknowledge this financial support of the SNSF. Finally, this work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon, also supported by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process. This research has made use of computing facilities operated by CeSAM data center at LAM, Marseille, France.

Appendix A SHINE blind test

Appendix B SHINE IFS observations

References

All Tables

All Figures

	Fig. 1 Distribution of the planet candidates found and missed as a function of the (S/N)tot given a perfect recombination of the images. The true negatives, all grouped at S∕N = 0, are not shown.
In the text

	Fig. 2 K-Stacker (S/N)KS derivative infunction of the orbit number sorted by S/N.
In the text

	Fig. 3 K-Stacker (S/N)KS divided by (S/N)tot of a perfect recombination as a function of (S/N)tot.
In the text

Fig. 4 Example of orbital parameters coming from one true positive solution of the SHINE blind test (see Sect. 3.2). In each subplot, each point corresponds to the parameters of one orbit found by K-Stacker with the color of its (S/N)KS value. In the diagonal, the histograms for each orbital parameter are shown. The dark cross corresponds to the mean of the K-Stacker orbits. The size of the dark cross is 1σ of each orbital parameter. The star is at the position of the orbital parameters used to inject the fake planet.

In the text

Fig. 5 Results of the K-Stacker algorithm when including an error on the stellar mass used by the algorithm to calculate the orbital positions at each epoch. The dots (resp. dashes) show the mean (resp. maximum)distance between the true position and the position found by K-Stacker in all the images as a function of the error on the mass. These dots are color coded to indicate the corresponding (S/N)KS. The black dashed line indicates the size of the FWHM of the instrumental PSF (i.e., K-Stacker has missed the planet at least in one image when a small line is above the dashed line).

In the text

	Fig. 6 Difference between the orbital parameters found by K-Stacker and the real values of injection, as a function of the error on the stellar mass.
In the text

	Fig. 7 Best recombined image resulting from the K-Stacker run on β Pictoris. At each epoch, the images are rotated and shifted to put the planet on its periastron position found by K-Stacker, and the frames are co-added. The planet b is detected at a (S/N)KS level of 24.5.
In the text

Fig. 8 Histograms and 2D diagrams of the β Pictoris b orbital parameters found by K-Stacker. Left, top, and bottom: scale of the orbital parameters. Right: scale of the histograms. 89 points in each 2D diagram corresponding to the 89 orbits with the highest (S/N)KS found by K-Stacker. The color of each point gives the (S/N)KS indicated at right. The dark cross indicates the mean value of the orbital parameters with their error bars. The red cross shows the higher probability density found by an MCMC technique in Lagrange et al. (2019a) converted in the K-Stacker referential. The origin of the t0 K-Stacker date is the 01∕12∕2014.

In the text

	Fig. 9 Best recombined image resulting from the K-Stacker run on HD 95086. At each epoch, the images are rotated and shifted to put the planet on its periastron position found by K-Stacker, and the frames are co-added. The planet b is detected at a (S/N)KS level of 9.97.
In the text

Fig. 10 Histograms and 2D diagrams of the HD 95086 b orbital parameters found by K-Stacker. Left, top, and bottom: scale of the orbital parameters. Right: scale of the histograms. The 129 points in each 2D diagram correspond to the 129 orbits of the higher (S/N)KS found by K-Stacker. The color of each point gives the (S/N)KS indicated at right. The dark cross indicates the mean value of the orbital parameters with their error bars. The red cross shows the higher probability density found by an MCMC technique in Chauvin et al. (2018) converted in the K-Stacker referential. The origin of the t0 K-Stacker date is 05∕05∕2015.

In the text

Fig. 11 Recombined images obtained with K-Stacker, when searching for additional companions around β Pictoris and HD 95086, two systems observed in the SHINE survey. At each epoch, the images are rotated and shifted to put the detection on its periastron position found by K-Stacker, and the frames are co-added. In each case, a bright spot can be seen near the coronagraphic mask (red arrow), which corresponds to the best solution found by the algorithm. In the case of β Pic, the corresponding (S/N)KS is 4.8. For HD 95086, the (S/N)KS is 4.4.

In the text