K-Stacker, an algorithm to hack the orbital parameters of planets hidden in high-contrast imaging. First applications to VLT SPHERE multi-epoch observations

Recent high-contrast imaging surveys, looking for planets in young, nearby systems showed evidence of a small number of giant planets at relatively large separation beyond typically 20 au where those surveys are the most sensitive. Access to smaller physical separations between 5 and 20 au is the next step for future planet imagers on 10 m telescopes and ELTs in order to bridge the gap with indirect techniques (radial velocity, transit, astrometry with Gaia). In that context, we recently proposed a new algorithm, Keplerian-Stacker, combining multiple observations acquired at different epochs and taking into account the orbital motion of a potential planet present in the images to boost the ultimate detection limit. We showed that this algorithm is able to find planets in time series of simulated images of SPHERE even when a planet remains undetected at one epoch. Here, we validate the K-Stacker algorithm performances on real SPHERE datasets, to demonstrate its resilience to instrumental speckles and the gain offered in terms of true detection. This will motivate future dedicated multi-epoch observation campaigns in high-contrast imaging to search for planets in emitted and reflected light. Results. We show that K-Stacker achieves high success rate when the SNR of the planet in the stacked image reaches 7. The improvement of the SNR ratio goes as the square root of the total exposure time. During the blind test and the redetection of HD 95086 b, and betaPic b, we highlight the ability of K-Stacker to find orbital solutions consistent with the ones derived by the state of the art MCMC orbital fitting techniques, confirming that in addition to the detection gain, K-Stacker offers the opportunity to characterize the most probable orbital solutions of the exoplanets recovered at low signal to noise.


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 the planet light that is drowned in the much brighter diffracted light from its host star. Jupiter and Earth like planets are about 10 8 to 10 10 fainter than their parent star in the visible Based on observations collected at the European Southern Observatory under programmes: 095.C-0298, 096.C-0241, 097. C-0865, 198.C-0209, 099.C-0127 band. However, thanks 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). They led to the discovery of the first planetary mass companions in the early 2000's, at large distances ≥ 100 au. The implementation of differential techniques, starting in 2005, en-Article number, page 1 of 13 arXiv:2004.12878v1 [astro-ph.EP] 27 Apr 2020 A&A proofs: manuscript no. LeCoroller abled the breakthrough discoveries of closer and lighter planetary mass companions like HR 8799 bcde (10,10,10 and 7 M Jup at respectively 14,24,38 and 68 au; Marois et al. 2008Marois et al. , 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 (GPI, Macintosh et al. 2014 ;SCExAO, Jovanovic et al. 2015;SPHERE, Beuzit et al. 2019) 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 −5 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 the past years has been the low occurrence rate of giant planets beyond 30 au (0.6 +0.7 −0.5 %, see Bowler 2016). Today, the GPIES and 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, submitted), suggesting that the bulk of the giant planet population is located between typically 1 and 10 au.
A prime goal of the future surveys will be also 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 like SPHERE+ (Boccaletti et al. 2020) aim at reaching contrasts of at least 10 −5 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 −6 at ≈ 60 mas) in 10 − 100 hours of exposure time, by providing an additional factor of 3 − 10 in contrast. In terms of observing strategy, for planets at less than 10 au from a star at 10 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 (higher contrast can be reached only 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 sub-optimal 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 PSF will be 4 to 5 times smaller than with the 10 m telescope class (VLT, Keck, Gemini, Subaru, LBT, LCO, etc.), and these ELTs will be used to search for planets at very small separations (below 10 au) already with the first-light instruments (Chauvin 2018a). In this case, the Keplerian motion could become non-negligeable in a matter of a few days only (Males et al. 2013).
Following similar principles previously applied to the search of new moons in solar systems like 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, and that 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 to 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 otherwise unreachable. This method can be used in addition to the 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 methods 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), by using simulated VLT-SPHERE observations, we have shown that when the total number n of available images is large enough to get √ n × (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 algorithms (PCA ADI and ASDI).
In Section 2, we describe the observations used in this paper, which come from the SPHERE SHINE survey (Chauvin et al. 2017a). In Section 3, we present the results of a blind test where K-Stacker was used to search for fake planets hidden in real 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 Section 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 Section 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 Section 6, a discussion on the strategy of future K-Stacker observations is done. Section 7 gives our final conclusions.

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 stars aimed 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. Even if the SHINE observations have not been organized for the K-Stacker 'philosophy' where we plan to observe less 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 to test K-Stacker for the first time in real conditions. SPHERE includes an extreme adaptive optics system, several types of coronagraphs, and three sub-systems (IRDIS, IFS and Zimpol). In this paper we use observations coming from the Integral Field Spectrograph (IFS) and the Infra-Red Dual-band Imager and Spectrograph (IRDIS), that were designed to cover the near-infrared range for an efficient search of young planets (Beuzit et al. 2019). In the SHINE survey, two 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 respectively optimized for the IRDIFS and IRDIFS-EXT modes (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 Section 3 were acquired using the APO1-ALC2 coronagraph. In Section 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) in order to constrain the orbital parameters of the known planets b around these targets.

Set-up 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 one 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 10 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 (seeing, AO performances, etc.). We create a total of 5 such fake K-Stacker runs, using a total of 50 different images from the survey. Table 1. Parameters used to inject the planet in the 30 K-Stacker runs of 10 observations of our blind experiment. We chose the same ranges than in paper I. Fake planets are 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 4 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 −6 , 4 × 10 −7 ]. 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 is repeated 6 times, in order 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 ratio of the injected fake planets. Note that, the S/N definition when using K-Stacker can be confusing, as we need to distinguish 3 different values: the S/N ratio 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 three following Groups: -Group I, composed of 13 runs for which (S/N) tot < 2 for the injected fake planets (or non injected fake planets). These planets can be considered as undetectable i.e. equivalent to no planet injected. -Group II, composed of 9 runs for which 2 < (S/N) tot < 12 for the injected fake planets (low-S/N regime) -Group III, composed of 8 runs for which (S/N) tot > 12 for the injected fake planets and therefore easily detectable. Following the principles of Paper I, after running the K-Stacker algorithm, we asked to an independent observer (not aware of the presence or not of fake planets in the images, of the injection phase, the S/N or the orbital parameters) to check the final K-stacker recombined solutions and to assign for each run one of the three flags: "no detection", "planet candidate", and "possible candidate, more observations required".

Blind check
Among the thirteen sets of Group I (i.e. (S/N) tot < 2), eight have been correctly flagged by the observer as no detection i.e. true negatives; The five other solutions were flagged as "possible Article number, page 3 of 13 Evolution of the SNR derivative on the best 100 orbits candidates, more observations needed". Due to the nature of the blind test performed, in which the same limited number of highcontrast observations were used multiple times, all the runs are not fully independent, and these 5 cases actually correspond to two truly independent candidates emerging from the noise. The typical K-Stacker S/N ratio 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) tot < 12), five solutions were flagged as true detections by the observer and were indeed true positives. Two solutions at (S/N) tot = 3.3 and (S/N) tot = 3.6 have been 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 than described above. Two other solutions at (S/N) tot = 3 and (S/N) tot = 4 have been 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 by the true positions of the fake planet only in 3 images (epochs) and not along the other last 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 the eight runs of Group III (i.e. (S/N) tot > 12) have been flagged by 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) tot > 12 are not shown. All these solutions with (S/N) 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) 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 Section 5). However, in every cases, 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.

Extracted orbital parameters
The result of the K-Stacker algorithm is a list of orbits sorted by signal to noise (see Paper I for detail). The signal to noise of the ≈ 50 − 100 first orbits is a relatively flat function before decreasing by steps (see example in Fig. 2). For our study of the orbital parameters, we keep only the orbits before this drop in (S/N) KS i.e. 68 first orbits in the example of Fig. 2. We have checked that for the true positives of the SHINE blind test the first K-Stacker solutions before the decrease of (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 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 3 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 Tab. A.1). Fig. 3 shows also that the K-Stacker signal to noise (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 well recovered the orbital parameters although the planet has travelled over a maximum of 15 − 38% of its total orbital period (see Fig. 4 and TableA.1).

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.

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 semi major-axis a, the orbital velocity is given by: Thus, a small error δM on the mass M of the star directly translates to an error δV orb on the velocity given by: This error on the orbital velocity will lead to a K-Stacker estimate of the position in each image which "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 : For a star at a distance d, the associated angular error is then: With 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 Figure 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 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 Equation 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, Figure 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 Figure 5. The reason for the existence of such a secondary maximum is not clear, and the grid sampling could play an important role here. 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 re-work 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 cur-rent paper but could be considered, if necessary, for example in the case of low mass stars where the mass is not always well known.

Tolerance on the true north error
The True North (TN) gives the absolute rotational orientation of observations. An error on the true north 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 , 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. , 2015. To check the consequences of a TN error on the K-Stacker algorithm, we simulated several observations taking into account the real distribution of the TN errors measured on SPHERE (see Tab. 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 FWHM of the PSF at the edge of the AO corrected area. At that level, this error should not affect the performance of K-Stacker in any noticeable way. And indeed, we found that K-Stacker had no problem to recover the planets, and gave the same final S/N ratios and orbital parameters, even when multiplying the SPHERE TN errors by a factor of 10.

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 DR2 for bright sources (mag < 14) is δ < 0.1 mas (Luri et al. 2018). For a star at 10 pc Error on a (a.u.) 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.
its central star, and then projects this position on the detector. The projected position is directly proportional to d −1 . Given that the AO corrected field 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. 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 moving by only a few PSFs over the 5 years of IRDIS monitoring.

β 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 ratio 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 was 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 sub-plots 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 which 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 0 found by Lagrange et al. (2019a) that is near the apoastron of the K-Stacker mean solution (see red and dark crosses for t 0 in Fig.  8). This ambiguity in the periastron/apoastron is reinforced by the small eccentricity and the incomplete coverage of the orbit. It will be solved with further monitoring.

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 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 signal to noise was improved by a factor 2.72 by K-Stacker. For a series of 8 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).

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 the known 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. A chromatic study shows that this bright speckle is a false positive (Desgrange et al. 2020, In preparation).
A&A proofs: manuscript no. LeCoroller   In the case of β Pictoris, the compact feature is found at (S/N) KS = 4.9 (Fig. 11) with the orbital parameters given at Table 5.
The semi-major 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 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 where the orbits found by K-Stacker passed well by the planet positions only at 3-4 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%.

Discussion on the 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 realised 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 characterising the physical properties of the planet (atmosphere composition, temperature, surface gravity, clouds detections, etc.); see for example the characterisation of the planets around the stars HR8799, β Pic, HD95086, 51Eri, HIP65426, PDS70 , Chilcote et al. 2017, De Rosa et al. 2016, Chauvin et al. 2017b. But, further work is required to develop a K-Stacker exoplanet statistical analysis tool : using a similar approach than with the multi-purpose exoplanet simulation system (MESS) algorithm (Bonavita et al. 2012), we should be able to inject numerous fake planets in series of observations in order 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 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 will also depend of 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 eleven exposures, a minimum of 11 * [(7/4.8) 2 − 1] ≈ 12 new observations are needed (perhaps a little less by taking care to observe β Pictoris c when it is fully out of the coronagraphic mask using radial velocity constraints) to reach (S/N) KS ≈ 7. Fig. 3 and Table A.1 shows that K-Stacker is able to recombine in a close to optimal way, 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 (i.e. in each observation, a parallactic rotation of at least one PSF FWHM at the position of the planet must be reached to subtract efficiently the instrumental speckles with ADI). The trade-off comes from the difficulty in setting up the observation strategy, the telescope overheads, and the 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 Article number, page 9 of 13 A&A proofs: manuscript no. LeCoroller

Beta pic
HD 95086 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.9. For HD 95086, the (S/N) KS is 4.4. split the observations over at least 10 − 20% of the orbital period of the searched planets to get accurate orbital parameters.

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 concluded 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 is 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 sub-stellar companions around these two stars. We found two bright features, corresponding to two possible planets orbiting within the orbit of the known companions. However, these two features are close to the coronagraphic mask, and we will show in a dedicated analysis that the one found around HD 95086 is a recombined "superspeckle". 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, θ 0 ) 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.  Table A.1. Orbital parameters of the injected and found planets (true positives) during the blind test described at Sect. 3. ∆ (pixels) is the mean distance between the injected and found positions of the planet in the images. S /N is the signal to noise computed using the orbits of injection (S/N) tot and the mean orbit found (S/N) KS . %T is the percentage of the covered orbit.