© The Authors 2023

1 Introduction

From an observational point of view, stellar activity poses a challenge to detections of planet around stars different from the Sun and to precise measurements of their parameters, especially with regard to the radial velocity and the transit photometry techniques. In stellar light curves (LCs), the out-of-transit flux is not typically flat but, instead, it is characterized by amplitude variations that are due to the intrinsic stellar activity, which is particularly enhanced in young stars (i.e., <1 Gyr). The latter do indeed show flux variations on many different timescales, with amplitudes up to 700% in a temporal range of 1–10 days (Cody et al. 2017; Cody & Hillenbrand 2018). These fluctuations in brightness are caused by several phenomena occurring on the stellar surface, primarily spots and faculae, that is, colder or brighter regions on the stellar surface, which can remove or add flux, depending on their location when the planet is transiting, thus affecting both the transit depth and shape, as well as the overall LC. Another source of noise for transit parameter determination are stellar oscillations, granulation, and flares, which contribute to increase the flux variability, consequently affecting planet detection. The effect of stellar activity on transit is more remarkable the smaller the planet is. In particular, young stars are characterized by fast rotation and complex brightness modulation, which makes them challenging targets for exoplanet research.

Despite these difficulties, observing young stars is important for multiple reasons: these young exoplanets can provide hints about the first stages of planet formation and evolution, helping us to put some constraints on current formation models. Understanding the mechanisms occurring in the early evolutionary phase of planet formation – such as migration and ionizing radiation from the host star (Baraffe et al. 2003), as well dynamical interactions with the protoplanetary disk, the other bodies in the system (e.g., the Kozai–Lidov mechanism; Shevchenko 2017), or with nearby stars (Ida & Lin 2010) – is fundamental in order to better interpret the final distribution of extrasolar planet population. Indeed, all these processes contribute to change both orbital and planetary parameters, such as mass, eccentricity, and distance from the host star, especially in the first phase of planet formation.

The aforementioned stellar activity in young stars can be of roughly on the same order of magnitude as that of the transit signal, in terms of both period and flux variation (Armstrong et al. 2015), thus masking the true planetary signals or producing false positives (Rodenbeck et al. 2018). This is the reason why most of the transit surveys are focused on older, photometrically quieter stars, exhibiting only small variations over long time scales (Rowe et al. 2014). Nevertheless, when looking for Earth-sized planets around solar-like stars, even this low variability of the stellar flux over long time scales can become problematic for the transit detection. Therefore, many filtering algorithms have been developed in order to correct for stellar activity in photometric LCs.

According to Hippke et al. (2019), in planet injection-retrieval experiments, these methods perform very well in quiet stars, reaching a detection efficiency of nearly 100%. However, this percentage decreases down to ~43.8% when considering a sample of young, extremely active stars observed with the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), meaning that less than half of the injected transit signals are recovered after detrending the highly variable LCs produced by these young stars. In their tests, Hippke et al. (2019) considered planets with a radius equal to half the Jupiter radius, hence, the recovery fraction of transit signals is expected to decrease even more when considering small planets such as mini-Neptunes and Super-Earths. This proves how challenging is searching for young exoplanets around young active stars.

Several algorithms have been proposed to specifically filter activity in the LCs of young stars in order to search for planetary transit signals, leading to the discovery of young exoplanets (e.g., Mann et al. 2016; Newton et al. 2019; Benatti et al. 2019; Nardiello et al. 2019; Thao et al. 2020). However, many more young exoplanets are needed for building reliable statistics.

The main purpose of this work is to identify and compare algorithms specifically developed for the filtering of stellar activity in young stellar LCs, in order to test their effectiveness prior to the search of transit events. We performed tests on realistically simulated LCs of both quiet and active solar-like stars, as these are the types that will be observed by the forthcoming ESA mission Planetary Transits and Oscillations of stars (PLATO; Rauer et al. 2014). In particular, we selected four custom-built algorithms: Notch and LOCoR¹ (N&L; Rizzuto et al. 2017), Young Stars Detrending² (YSD; Battley et al. 2020), K2 Systematics Correction³ (K2SC; Aigrain et al. 2016) and VARLET and PHALET (Grziwa & Pätzold 2016).

These have been recently successfully applied by independent groups on Kepler (Borucki et al. 2010), K2 (Howell et al. 2014) or TESS LCs, specifically for the search of planets around young, active stars. Along with the above-mentioned algorithms, we also included in the comparison two of the algorithms selected as the best performing ones for less active stars according to Hippke et al. (2019), namely, “Tukey’s biweight” and “Huber spline,” as implemented in the Wōtan⁴ Python package. This paper is organized as follows. In Sect. 2, the PLATO mission is introduced, along with the description of the samples of PLATO simulated LCs and the criteria for transit detection. A brief summary of the algorithms tested by Hippke et al. (2019; i.e., general-purpose algorithms) and the analysis of their performance on a series of different transit signals injected into PLATO LCs is presented in Sect. 3. In Sect. 4, the custom-built algorithms mentioned above are introduced. Then, in Sect. 5, several injection-retrieval tests are performed on both individual quarters (Sect. 5.1) and entire LCs (Sect. 5.3). Moreover, we investigate some peculiar cases: eccentric hot Jupiters (Sect. 5.4.1) and planets with an orbital period equal to the rotation period of their host star (Sect. 5.4.2). Finally, in Sect. 6, we discuss the main results and future perspectives.

2 PLATO simulated light curves and transit detection

PLATO⁵ (Rauer et al. 2014) is the third medium-class mission in ESA’s Cosmic Vision program. Scheduled for launch in 2026, PLATO has the main aim of detecting and characterizing, in terms of density and age, Earth-like planets around solar-like stars up to the habitable zone. PLATO will carry out high-precision photometric observations on a large sample of stars in order to characterize both the exoplanets and their host stars, respectively, through transit detection and asteroseismology. Since PLATO is designed to observe a broad field of view (FoV⁶; Nascimbeni et al. 2022), similarly to the NASA all-sky survey mission TESS, a multi-telescope approach was adopted. Indeed, PLATO is equipped with 24 “normal” cameras or N-CAM (i.e., 25 s readout cadence), organized into four groups of six cameras, each with its own CCD-based focal plane array, and two “fast” cameras or F-CAM (i.e., 2.5 s readout cadence) dedicated to the observation of bright stars through the use of two bandpass filters (one per camera) acting as fine guidance sensors. All the cameras are refracting telescopes with an aperture diameter of about 120 mm, pointing at different parts of the sky (with some overlapping areas), covering a total FoV of more than 2200 deg² (Pertenais et al. 2021; Magrin et al. 2020). Each camera is equipped with four charge-coupled devices (i.e., 104 CCDs in total) with 4510 × 4510 pixels of size 18 μm (edge length). In preparation of the future PLATO space mission, several research groups implemented simulated LCs as those expected to be obtained by the PLATO telescope. According to the current observation strategy, PLATO will point the same portion of the sky for at least 2 yr, but slightly adjusting the spacecraft position every 3 months (i.e., every quarter), rotating 90° around the line-of-sight. Hence, considering the total time-series datasets, between one quarter and the next, a jump in the LC is expected, due to the fact that the pixel position on the camera where a given target is observed can change between a quarter and another.

2.1 Simulated light curves

For our analysis, we employed two samples of 100 simulated LCs provided by the Light curve Stitching Working Group (LSWG), which is investigating the best method to model the flux variations between quarters of PLATO LCs. Each LC encompasses eight quarters with a duration of ~88 days and gaps between quarters of 1–3 days, covering 2 yr in total. The LCs were originally produced at a sampling of 25 s. Since such a high cadence was not needed for the transit detection, in order to save computational time and keep the file size manageable, the curves were binned to 600 s cadence.

The LSWG produced 100 LCs simulating PLATO-like data of bright stars (8 ≤ V ≤ 11 mag)⁷ with a photometric amplitude variability in the range 0.08–1.71%, representative of PLATOlike systematics, as obtained by the PLATO Solar-like Light curve Simulator (PSLS) 1.2⁸ (Samadi et al. 2019). We note that for this particular set of simulations we did not use any realistic pointing direction, that is, the noise values are representative of what PLATO will observe, but there is no correlation with the distribution of stars in the sky.

As the mission is still under development, many technical details such as the noise levels, the duration of the runs (in this work: 88 days), and the duration of the gaps between the quarters (in this work: between 1–3 days) are based on requirements or realistic assumptions, while real in-flight performance might be different. The instrumental systematics and random noise for the noise-only LCs were simulated assuming point spread function (PSF) photometry (as will be done on the ground for the bright stars in the P1 sample). The final parameter that controls the noise in PSLS is the stellar magnitude, which is reported in Tables B.1 and B.2, along with the full set of input parameters for the LCs used in this work.

For the activity signals, we required a physically motivated model where individual parameters such as the rotation rate and the spot coverage could be controlled by the user. This is only partially accomplished by the stellar activity component in PSLS, which consists of a Lorentzian in Fourier space with no explicit periodicity (see Sect. 5.2 of Samadi et al. 2019 for more details). Hence, we switched off the variability components of PSLS and manually inserted activity signals and oscillations afterwards. Therefore, onto these noise-only LCs, we superimposed the activity signals of solar-like “quiet” stars with P_rot ~ 1–50 days extracted from LCs simulated in Aigrain et al. (2015), which included Sun-like butterfly patterns, activity cycles, spot evolution (i.e., emergence and decay)⁹, and differential rotation, on a temporal range of 1000 days and representative of the Kepler data of solar-like stars. We refer to this dataset as the “quiet sample.”

The “active sample” required to test the different algorithms was based again on the same LCs produced by the LSWG, this time with the addition of a strong quasi-periodic activity signal with a rotation period shorter than 10 days¹⁰, again extracted from the simulations performed by Aigrain et al. (2015). Two examples, representing the LCs of a quiet and an active star respectively, as simulated by the LSWG are shown in Fig. 1 (top) and (bottom), respectively.

Finally, we did not attempt to remove the instrumental long-term systematics which are very obvious in the top panel of Fig. 1, as the algorithms that will address those residuals in the PLATO pipeline are still under development. The LCs employed in this analysis should not be considered as representative of PLATO final data products, which indeed will be corrected for instrumental systematics¹¹ but should rather be considered a worst-case scenario in which intermediate data products are employed for the analysis.

Fig. 1 PLATO simulated LC #0 of the sample of quiet stars (top) and active stars (bottom). Different colors represent different quarters. The systematics in the simulated LCs represent a worst-case scenario, as their correction will be addressed in the final PLATO data products.

2.2 Transit modeling

Transits were modeled using the batman¹² Python package (Kreidberg 2015), which is a Python implementation of the transit model from Mandel & Agol (2002). Since we wanted to measure the performances of the algorithms on different populations of exoplanets, the orbital parameters and planetary properties were varied from case to case, and they are described along the paper accordingly. The injection was performed by simply multiplying the stellar LC by the transit model (normalized to unity outside the transit), that is, spot-crossing events and other effects altering the shape of the transit were not included. The model did not include other planet-related effects other than the transit, for instance, phase curve variations, secondary eclipses, Doppler boosting, and ellipsoidal variations were not included. After the transit injection, these LCs were detrended with selected filtering algorithms (see Sects. 3 and 4). The transit detection on the resulting filtered flux was then performed by the Transit Least Squares (TLS) algorithm (Hippke & Heller 2019), which has shown a higher detection efficiency compared to other transit detection algorithms, such as Box-fitting Least Squares (BLS; Kovács et al. 2002). The main improvement of TLS with respect to BLS is that it employs a real planetary transit shaped model when searching for transit signals in the phase-folded corrected LC and not simply a less realistic boxlike shape. The TLS default template of the transit curve used for transit detection in this work, assumes circular orbits (i.e., a null eccentricity), with an impact parameter of 0 (i.e., an inclination of 90 degrees). Unfortunately, we realized at a later stage that the filtering algorithms may modify the transit signal to a point that there is no real advantage in using a more accurate shape for the transit; thus, in the future, we advise to use more than one detection technique in order to maximize the scientific return.

TLS computes the Signal Detection Efficiency (SDE) as the detection statistics for finding the strongest signal in the periodogram consistent with a planetary transit. The SDE distribution as a function of the orbital period, is obtained in the same way as in Kovács et al. (2002): (1)

where SR(P_orb) is the signal residual between the model and the data, computed from the distribution of minimum χ⁽²⁾. In particular: 〈SR(P_orb)〉 is the arithmetic mean, σ(SR(P_orb)) is the standard deviation, and SR_peak is the peak value, which (by definition) is equal to 1. Once that the SDE distribution has been computed, the period corresponding to the highest SDE value is assumed as the orbital period of the transit signal. In our analysis, the planet was considered to have been recovered if the detected period matched the injected one within 1%. Moreover, for a more realistic transit search, a detection SDE threshold of 7 was applied, as done by Hippke & Heller (2019). This means that all signals with SDE below this threshold were not considered as correctly identified, even if the recovered orbital period matched the injected one within 1%. We also made use of the SDE of the recovered period as a metric for the comparison of different algorithms. We are aware that the SDE threshold should depend on the properties of the light curves and that a single value may not be appropriate (see for example Grziwa et al. 2012), but we believe it still represents a reasonable choice given the goals of this work. We stress again the fact that our choice is not representative of the final implementation of the PLATO pipeline, and it is only valid in the framework of this specific work. Finally, we note that we did not produce any comparison sample, that is, a set of LCs without injected planets, as the goal of this work is to test the ability of different algorithms in recovering bona fide planets rather than determining the PLATO’s planet yield.

In the next section, the most successful general-purpose filtering algorithms are briefly described and applied to the samples of PLATO simulated quiet and active stars.

3 General-purpose algorithms

Hippke et al. (2019, hereafter H19), performed planetary injection-retrieval tests on a sample of 316 highly variable stars observed with TESS, obtaining a global detection efficiency lower than 50% when applying algorithms that were rather successful in detecting small planets around quiet stars. They injected transit signals with the following properties: a planetary radius equal to half the Jupiter radius (R_p = 0.5 R_J), with period randomly drawn from a uniform distribution in the range 1–15 days, and an impact parameter and eccentricity both fixed to zero. We reproduced the same experiment on the samples of both quiet and young active stars of the PLATO simulated LCs, in order to test the performances of the general-purpose filtering algorithms used by H19 (already implemented in the Wōtan python package) on PLATO LCs. In particular, the main aim of this analysis is to investigate how the cadence and the duration of the observations of PLATO can influence the efficiency of the algorithms presented by H19 or if these algorithms are not adequate for the detection of young exoplanets by design. Specifically, H19 tested the following algorithms¹³: sliding median, Tukey’s biweight, Huber spline and Lowess.

In the statistical literature, the sliding median (window-size of w = 0.7 days) is classified as a “scatterplot smoother” that uses a rectangular (i.e., symmetric and with uniform weights) window kernel of length, w, centered on the time, t(x_i), for each data-point, x_i.

The Tukey’s biweight (w = 1.0 or 0.25 days for the quiet and active sample, respectively) is a robust location estimator similar to the ordinary least-squares method except that the weights are not constant but depend on the distance from the midpoint, while the loss function¹⁴ is given by: (2)

where a are the esiduals (i.e., the difference between the observed data and the midpoint) and c is a constant called tuning parameter, usually with a value of ~5 (default in Wōtan).

The Huber spline (w = 0.3 days) is a generalization of the maximum likelihood estimation, which is characterized by the following approximated loss function: (3)

where c = 1.5 as the Wōtan default option.

Finally, Lowess regression (w = 1.0 day) consists in a local polynomial regression method which works by fitting a low-order polynomial to a subset of the data (defined by a certain window-size) at each point along the x-axis using weighted least-squares regression. In this way, more weight is given to the points closer to the data point being estimated.

We will refer to these algorithms as general-purpose as they were not specifically developed to deal with prominent stellar activity. In order to test the efficiency of the aforementioned general-purpose algorithms on PLATO LCs, we performed planetary injection-retrieval tests that were as similar as possible to those done by H19. For this analysis, we treated each quarter as an independent LC, that is, the test was performed on 800 LCs with 88 days of duration. Planetary transits with the characteristics described above were synthesized with batman and then injected into the single quarters of both the quiet and the active LCs of PLATO. We injected a different planet into each quarter, so that in total, 800 transit signals were analyzed. For the stellar density, ρ_*, and for the limb darkening coefficients, u₁ and u₂, we adopted the values reported in the tables describing the simulated LCs, provided by the LSWG (shown in Appendix B). Afterward, the LCs with injected transits were filtered with the Wōtan algorithms with the same window-sizes adopted by H19. Finally, we analyzed the filtered fluxes by searching for planetary transits with TLS.

The general-purpose algorithms achieved unexpectedly excellent results on both samples of the PLATO simulated LCs, when compared to the results presented in H19. The percentage of recovered transits with TLS after the filtering with these algorithms is reported in Fig. 2, where each row represents the results of a different algorithm, and each column refers to a different sample of LCs. Specifically, in the first column the results obtained by H19 on TESS active stars are reported, whereas in the second and third columns, the results from the PLATO quiet and active sample are shown, respectively. In the last row, the recovery percentage obtained combining all methods is also shown. As can be seen from this figure, the recovery efficiency of the tested general-purpose algorithms on both samples of PLATO stellar LCs is extremely high, recovering 100% of the injected signals in almost all the cases. The only exception is the biweight algorithm with a window-size of 0.25 days, which recovers 98.88% and 98.62% of the transit signals in the quiet and active sample, respectively. In a reanalysis of the LCs of both samples with the biweight algorithm and a window-size of 1 day – a good approximation of the best value proposed by H19 (three times the duration of the transit) for orbital periods up to 15 days – the TLS recovers the totality of the injected transit signals.

An example of a planet missed by TLS after filtering the LC with the biweight with w = 0.25 days is given in Fig. 3, where the raw and the detrended flux of quarter #3 of LC #9 of the active sample are shown. In this quarter, a planet with an orbital period (P_orb) of 12.60 days and a scaled planetary radius (R_p/R_*) of about 0.056 is injected. In each panel of this figure, the stellar variability model obtained by means of different algorithms is highlighted in different colors: red for the median, green for the Huber spline, orange for the Lowess, blue for the biweight (specifically, dark blue for a window-size of w = 0.25 day and light blue for w = 1.0 day). It can be easily observed that the biweight with the shortest window fails in correctly filtering the LC, so that in the corresponding corrected flux, almost all the transits disappear. Nevertheless, by increasing the window-size up to 1.0 day, transit features are perfectly recognizable after the filtering and the transit is correctly detected.

Considering the SDE of the recovered transits computed by TLS in the active sample, in the boxplot in Fig. 4, the main properties of the SDE distributions are summarized. It can be immediately visualized that the SDE for most of the algorithms ranges between ~30 and ~76, whereas the biweight with the shortest window-size reaches much lower values. Similar results are obtained from the analysis of the quiet sample of PLATO simulated LCs.

Although we cannot exclude that the level of activity in the PLATO LCs is lower than those observed in the TESS dataset analyzed by H19, our analysis shows that the low efficiency general-purpose algorithms for the filtering of active stars derived by H19 are likely driven by the characteristics of the LCs, that is, the longer time span and the noise properties. For this reason, we decided to include the two¹⁵ best-performing algorithms in our comparison and in the analysis presented in H19. Overall, the biweight with w = 1.0 day and the Huber spline with w = 0.3 day appear to be the algorithms with the highest values of SDE and the lowest standard deviation, with regard to the quiet and active sample, respectively. Therefore, these two algorithms are tested also on the other injection-retrieval experiments performed, considering also smaller planets, as explained in the next sections.

Fig. 2 Recovery fraction of injected transit signals similar to those of Hippke et al. (2019) injected into simulated PLATO stellar LCs of quiet (second column) and active (third column) stars, compared with the recovery fraction of similar signals injected into a sample of active stars observed by TESS (first column) and analyzed by H19. The biweight with a window of 1.0 day was not tested by H19 on TESS active LCs.

4 Custom-built algorithms

As previously mentioned in Sect. 1, four independent research groups have recently developed algorithms specifically for the filtering of activity in the LCs of young stars observed by K2 and TESS. In this work, they were implemented in Python and applied over the two samples of PLATO simulated LCs. A brief description of the selected algorithms is given below.

N&L, developed by Rizzuto et al. (2017), is a two-step pipeline, called Notch and LOCoR, respectively. Notch works by fitting two models over a moving data window: a simple quadratic polynomial and a polynomial with a box-shape transit notch, in order to identify the possible presence of a transit. Then, the model with the lowest Bayesian Information Criterion (BIC) is adopted. The selection of the window-size employed for the filtering of PLATO LCs is summarized in Eq. (4) and it depends on the stellar rotation period P_rot: (4)

For P_rot ≤ 2 days, Notch is unable to completely remove rotational systematics, therefore the LOCoR algorithm should also be applied. The latter is based on the LOCI algorithm of Lafrenière et al. (2007) and it consists of modeling each individual rotation as a linear combination of other rotations of the same star in the dataset.

YSD, developed by Battley et al. (2020), is a Python pipeline which employs Lowess smoothing regression to model stellar variability, using the standard tricube as a weighting function: (5)

where the weight w at each data point, x, is given by the distance, d, from the point on the curve being fitted, scaled to lie in the range [0,1].

Therefore, the detrending pipeline is carried out by cutting the outliers¹⁶ (peaks and troughs) first, and then by estimating the variability trend using the Lowess function implemented in the Python module statsmodel (Seabold & Perktold 2010) with frac = 0.012 when considering one quarter only, meaning that a fraction of 1.2% of the data has been used to estimate each value of the variability model. When filtering the entire PLATO LCs, the frac parameter is set to 0.0015 instead, because it results in better performances.

The outliers are identified with the function find_peaks in the scipy.signal (Virtanen et al. 2020) package¹⁷, using a prominence of 0.001 and a width of 20 data points. These parameters should be adjusted for more complex LCs. For more information about the outlier cutting see Sect. 2.5.2 of Battley et al. (2020).

K2SC, developed by Aigrain et al. (2016), is one of the filtering algorithms that makes use of Gaussian processes (GPs: Rasmussen & Williams 2006) to model both instrumental systematics and astrophysical variability. Overall, GPs are non-parametric methods which describe a dataset of N data points by evaluating correlations between them through a kernel or covariance function. The latter describes how each point is related to all the other points, and this relation is expressed through a N × N matrix, called a covariance matrix. To speed up computation time, we employed the covariance functions from the celerite2¹⁸ (Foreman-Mackey 2018) package in our implementation of Gaussian processes for the PLATO LCs.

Activity is modeled as a mix of two stochastically driven and damped harmonic oscillators (SHO), describing two modes in Fourier Space at the stellar rotation period, P_rot, and its first harmonic, respectively, and implemented in the code as the RotationTerm kernel. It was shown that this kernel is able to reproduce a large range of stellar variability in time-series dataset, from stellar rotation to pulsations (e.g., Haywood et al. 2014; Rajpaul et al. 2015; Uttley et al. 2005). Besides the RotationTerm (K_rot), describing the time-dependent variability of the light curve, the JitterTerm (K_wn) implemented in celerite2 is added to model the non-periodic variability or white noise, This second term is given by K_wn,ij = σ²δ_ij, where δ_ij is the Kronecker delta function. Therefore, the final covariance function, defining each element of the covariance matrix with indices (i, j), is given by: K_ij = K_rot,ij + K_wn,ij. The filtering has been done by means of the exoplanet¹⁹ package (Foreman-Mackey et al. 2021), fitting the data and imposing priors through the PyMC3²⁰ (Salvatier et al. 2016) statistical distributions, following Aigrain et al. (2016) as much as possible, even if the adopted kernels were quite different from the original pipeline of K2SC. For this reason, from now on, we will refer to this filtering algorithm as GPs and not K2SC. Despite the general name adopted in this work (GPs), it is important to highlight that Gaussian processes can be tested with many different kernels. In this work, we investigate only the specific kernel described above and a different kernel would likely return different results.

We adopted the following priors over the hyperparameters of the covariance function, with distributions defined in PyMC3: (1) the primary period of variability, fixed to the rotation period P_rot given in tables describing the simulated LCs, provided by the LSWG; (2) the standard deviation of the process σ, equivalent to the amplitude of the variability, with uniform prior in logarithmic space between 10⁻⁷ and 10⁰ in unit of normalized flux (i.e., flux divided by its median); (3) the quality factor of the secondary oscillation minus one half Q₀, with normal prior with mean 0.0 and standard deviation 2.0 in logarithmic space; (4) the difference between the quality factors of the first and second modes dQ, with normal prior with mean 0.0 and standard deviation 2.0 in logarithmic space; (5) the fractional amplitude of the secondary mode compared to the primary f, with uniform prior between 0.1 and 1.0.

VARLET, developed by Grziwa & Pätzold (2016), is a wavelet-based filtering pipeline divided into three main steps. Firstly, the raw flux is decomposed in wavelets at multiple levels through the discrete wavelet transform (DWT). Then, a thresholding function is applied to the details coefficients obtained from the previous multilevel decomposition. Specifically, a soft-thresholding procedure has been adopted (Donoho & Johnstone 1994), meaning that the corrected version (Y) of the original detail coefficients (X) is given by: Y = sign(X) * f(|X| − thr), where thr is the chosen threshold and f(|X| − thr) is defined by the following equation: (6)

where the threshold is computed through the penalized criterion provided by Birgé & Massart (2007). Finally, the denoized light curve is reconstructed through the inverse discrete wavelet transform (IDWT) using the corrected coefficients. Therefore, the reconstructed LC, which should represent stellar variability, is subtracted from the raw flux, thus obtaining a residual LC containing only the transit signals and the white noise.

Originally written in MATLAB, the pipeline has been reimplemented in Python by adopting the PyWavelets²¹ package (Lee et al. 2019). In the original pipeline, the VARLET filter is usually used with a standard BLS which does not take the transit depth or the shape (other than boxlike) into account. In this work, the search for transit signals in the filtered LC is instead carried out by means of the TLS.

All these algorithms have been tested together with the biweight and Huber spline algorithms described in Sect. 3.

Fig. 3 Raw flux (left panel) and detrended flux (right panel) of the 3rd quarter of LC #9 of the sample of active stars with an injected transit signal with Porb = 12.60 days and Rp/R*=0.056. Data-points are shown in black, whereas different colors represent the photometric variability model inferred from algorithm: Median (red), biweight with a window-size of 0.25 days (dark blue) and 1.0 day (light blue), Huber spline (green) and Lowess (orange).

Fig. 4 Boxplot of the signal detection efficiency for the general-purpose algorithms from Wōtan tested on the active sample with injected transit signals with Rp ~ 0.5 RJ and Porb ∈ [1–15] days, moving on circular orbits. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines the median. The black circles in the biweight (1.0 day) and the lowess (1.0 day) cases represent data that extend beyond the whiskers.

5 Injection-retrieval tests

In order to test the transit detection efficiency of the filtering algorithms described above, a series of planetary transit signals were synthesized and injected into the PLATO simulated LCs before filtering.

Transits representing planets with different scaled planetary radii (R_p/R_*) in the range 0.01–0.1, orbiting their host star in circular orbits (i.e., zero eccentricity), with an orbital period (P_orb) between 0.75 and 40 days and an impact parameter (b) drawn from a uniform distribution in the range [0–1], have been synthesized with the batman package and then injected into single quarters of both the quiet and active sample. We sampled the scaled planetary radius R_p/R_* from uniform distributions in the following ranges (𝒰(a, b)): 1) Jupiter-size: R_p/R_* from 𝒰(0.05, 0.1); 2) Neptune-size: R_p/R_* from 𝒰(0.03, 0.05); 3) Earth-size: R_p/R_* from 𝒰(0.01, 0.03).

Each planetary type (Jupiter-, Neptune-, and Earth-size) was analyzed independently from the other, that is: each LC was injected and analyzed three times (once per each planetary type), without mixing the results. Examples of the first quarter of LC #0 of both the quiet and active samples of PLATO simulated LCs, with injected transits for a Jupiter-, Neptune-, and Earth-like signal, are shown in Figs. 5–7, respectively.

As previously done, we retrieved the limb darkening coefficients and the stellar density (ρ_*), which is in the range 0.8–3.0 g cm⁻³, taken from tables provided by the LSWG. The scaled semi-major axis (a/R_*) were derived from the orbital period and the stellar density through Kepler’s third law.

5.1 Individual quarters

We performed the first transit injection-retrieval test considering each quarter of a simulated LC independently from the others, for a total of 800 tests for each planetary type (Jupiter-, Neptune-, and Earth-size type) as each of the 100 LCs comprises eight quarters. Although the eight quarters from the same LC are sharing the same stellar parameters, the different instrumental systematics and the quasi-periodic and non-stationary nature of stellar activity characterize each quarter in a different way, thus justifying our choice in the framework of testing the efficiency of algorithms. A given transit model is then injected in both the quiet and active samples, while keeping the samples separated. It follows that each efficiency rate is computed over 800 tests.

Each LC is then filtered using the algorithms described in Sect. 4 and the two general-purpose algorithms (i.e., the biweight and the Huber spline). The percentage of recovered transits with TLS after the filtering with these algorithms is reported in Figs. 8 and 9, for the quiet and active sample, respectively. Each row represents the results of a different algorithm, whereas each column refers to a different planetary type, as in the legend.

First of all, it can be noticed that in the quiet sample all the algorithms perform quite well for Jupiter- and Neptune-sized planets, as all filtering algorithms lead to transit recovery rates with TLS higher than 99.0%, except for VARLET which achieves 97.6 % for the Neptune class. For Earth-sized planets (i.e., R_p/R_* ≤ 0.03), all algorithms achieve retrieval rates higher than 90% except for VARLET (86%) and GPs (72.8%).

With regard to the active sample, a slightly decrease in the number of recovered transits has been observed, as expected, due to the much higher photometric variability of these LCs. Nevertheless, among the custom-built algorithms, N&L results to be the best-performing algorithm for all planetary types in both samples, recovering 96.87% and 90% of small planets in quiet and active stars, respectively. Among the general-purpose algorithms, the biweight and the Huber spline have performances extremely close to those of N&L in the quiet sample and the active sample respectively. Among the algorithms with a lower detection efficiency, we find the VARLET code for what regards the larger planets, that is, Jupiters and Neptunes, while the GPs slightly underperform with respect to the other methods for what regards the shallow transits generated by small planets.

H19 observed that GPs tend to absorb the transit signal by modeling it as part of the photometric variability of the light curve, thus reducing the transit depth as well as the sensitivity in the subsequent transit search in the filtered flux. This behavior was noticed also in the tests on PLATO LCs, not only with the GPs, but also with VARLET. An example is given in Fig. 10, where the Earth-like transit points in LC #22 are shown, with overplotted variability models inferred from different algorithms. It can be easily noticed that the filtering models determined by either the GPs (third top panel on the right, fuchsia line) or the VARLET model (last top panel on the right, violet line) exhibit a dip in correspondence of the transit, which is then shallower in the detrended flux (bottom panels), thus reducing the sensitivity and the SDE in the subsequent transit search. In the case of VARLET, this behavior is observed also for giant planets with large transit depths.

In general, we expect higher values of the SDE in the detection of planets with larger size (i.e., deeper transits) or shorter orbital periods (i.e., with a higher number of transits). This behavior was observed for all the filtering algorithms – except for VARLET, which generates a SDE distribution with similar low values regardless of the planetary size, as a consequence of systematically reducing the transit depths of all transits.

In particular, VARLET performs weakly in terms of SDE, achieving a maximum value of ~54, whereas all the other algorithms show similar performances, with SDE values up to 90 for the large planets. Similar results are obtained from the filtering done by the biweight and the Huber spline. Figures 11 and 12 show a performance summary of the analyzed methods on the quiet and active sample, respectively. With regard to the SDE of the general-purpose algorithms, it can be noticed that the biweight performs slightly better than the Huber spline in the case of quiet stars, whereas the situation is reversed in the case of active stars.

Finally, we note that the SDE of the detected transits on the filtered LCs of the active stars are similar to those of the quiet sample. This is an encouraging result, since it suggests that the activity filtering algorithms perform equally well for the two datasets, even for stars that are very different from the photometric point of view. This in turn suggests that PLATO will be able to detect a large number of planets with very different properties both in terms of size and periods, also around young or photometrically active stars.

Fig. 5 First quarter of the PLATO simulated LC #0 of the quiet (left) and active (right) sample, with an injected Jupiter-size planet. In black, we show the normalized flux without transits and in red, the transit model generated with batman.

Fig. 6 First quarter of the PLATO simulated LC #0 of the quiet (left) and active (right) sample, with an injected Neptune-size planet. In black, we show the normalized flux without transits and in blue, the transit model generated with batman.

Fig. 7 First quarter of the PLATO simulated LC #0 of the quiet (left) and active (right) sample, with an injected Earth-size planet. In black, we show the normalized flux without transits and in green, the transit model generated with batman. In the active case, the transits are hard to distinguish, therefore a zoom-in on the first transit is shown. The second transit occurs at about 70 days.

Fig. 8 Recovery fraction of injected transit signals in simulated PLATO stellar LCs of the sample of quiet stars.

Fig. 9 Recovery fraction of injected transit signals in simulated PLATO stellar LCs of a sample of active stars.

5.2 Caveats in transit detection

An in-depth analysis of the parameter space of the recovered planets shows that Jupiter and Neptune-size planets are mostly recovered regardless of their orbital period. On the other hand, Earths and super-Earths are harder to identify, especially at long orbital periods, since the number of transit signals decreases with increasing P_orb. This is to be expected since Earth-like planets produce shallow transit signals in the LCs, which can be overcome by the variability amplitude. An example of this trend is shown in Fig. 13, where the scaled planetary radius of injected planets is plotted as a function of the injected orbital period. Blue dots represent planets recovered by the TLS in the LCs of the active sample filtered by N&L (top) or Huber spline (bottom).

Furthermore, considering the orbital and planetary parameters of the recovered planets measured by TLS, it is found that filtering algorithms allow to correctly detect the mid-transit time and to measure a very precise estimate of the planetary orbital period with small uncertainties. However, this is not valid anymore when considering the scaled planetary radius R_p/R_*, which is often underestimated, especially for large planets (i.e., R_p/R_* ≥ 0.05). As an example, in Fig. 14, the orbital period and scaled planetary radius recovered with TLS after the filtering of the LCs of the quiet sample by N&L are compared to the corresponding injected values, with different colors representing different planetary size.

A possible explanation could reside in the TLS default template, which assumes circular orbits with a null impact parameter, as mentioned in Sect. 2. Whenever the transit shape strongly differs from the default, the recovered scaled planetary radius may be underestimated with respect to the original injected value, while also affecting all the other related parameters, such as the semi-major axis and the transit duration. This occurs, for example, in the case of grazing transits (i.e., high b), eccentric orbits, or when the filtering algorithm distorts the original transit, as in the case of GPs or VARLET. This is a known issue of the VARLET algorithm, which however does not represent an obstacle as the search for transit signals is carried out by means of the BLS algorithm, which does not take the depth of the fitted box into account.

Nevertheless, in cases where the transit shape differs from the TLS default template, the correct parameters can be retrieved by performing an additional analysis in order to fit the correct transit model to the LC. For example, in the case of grazing transits or eccentric orbits, a Markov chain Monte Carlo (MCMC) analysis on the filtered flux is sufficient to recover the orbital parameters correctly, as shown in Appendix A. When the orbital parameters are affected by a deformation of the transit shape, it is still possible to recover the original transit parameters by repeating the filter process after masking out the in-transit data points, which can be easily done after the transit has been detected. Therefore, we stress out the fact that the filtered LC used for the transit detection cannot be directly used for the planet characterization. An in-depth analysis is necessary for the precise measurement of the orbital and planetary parameters, either after reprocessing the LC while making use of the knowledge of the planet or by simultaneously fitting for stellar activity and planetary parameters.

Fig. 10 Raw flux (top) and detrended flux (bottom) of the first quarter of LC #22 of the quiet sample, zoomed in on the transit points of the injected Earth-size planet. Data points are shown in black, whereas different colors represent the photometric variability model inferred from different algorithms: N&L (orange), YSD (red), GPs (fuchsia), and VARLET (violet).

Fig. 11 Boxplot of the SDE for the algorithms listed on the y-axis, tested on the quiet sample with injected transit signals as described in Sect. 5.1. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines the median. The black circles in the VARLET case represent data that extend beyond the whiskers.

Fig. 12 Boxplot of the SDE for the algorithms listed on the y-axis, tested on the active sample with injected transit signals as described in Sect. 5.1. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines the median. The black circles in the VARLET case represent data that extend beyond the whiskers.

5.3 Entire light curves

One of the advantages of PLATO with respect to other space missions, such as TESS and Kepler, is the continuous observation of bright targets for at least two consecutive years, according to the current observation strategy. Such a long temporal baseline should allow to detect also the shallower transits of long-period small planets. In our previous analysis (i.e., Sects. 3, 5.1 and 5.2), we analyzed only one quarter of data per planetary system, so it is possible that some small planets were not detected because there were not enough transits to allow for their detection, regardless of the filtering algorithm employed. To test this scenario, we selected the Earth-size planets (i.e., R_p/R_* ≤ 0.03) left undiscovered after the N&L, Huber spline, and YSD filtering (i.e., the three algorithms with the highest recovery efficiency; see Fig. 9) and then injected their transits in the remaining seven quarters of the corresponding (active) LC. We then ran again the filtering algorithms and the TLS onto the entire LCs²². In this way, some of the previous missed transits are now detected whereas others are still missed, due to the extreme variability and high background noise of the active stars. Specifically, the percentage of recovered Earth-sized transits for N&L increases up to 94.025%, while for YSD, it reaches 96.55% and for Huber spline, it is 98.075%.

Surprisingly, the YSD and Huber spline algorithms seem to achieve better results than N&L when considering a much larger dataset²³, recovering even some of the small planets still missed by N&L.

An example is given in Fig. 15, where the raw, detrended, and phase-folded fluxes of LC #58 are shown in the upper, middle, and lower panels, respectively. The injected transit represents a small planet with a scaled planetary radius of R_p/R_* = 0.0146 and an orbital period of P_orb = 32.83 days. In looking at the filtered fluxes (middle panel), it is quite hard to distinguish transit signals (marked by the vertical blue lines), especially in the case where the YSD model has been applied, due to a lot of residual noise left after the filtering. Nevertheless, when computing the phase-folded flux (bottom panel), in the N&L case, the transit feature is barely distinguishable; in addition, it is so shallow that the TLS algorithm does not identify it as the strongest periodic signal, thus missing the planet. On the other hand, despite some residual variability, the transit feature is definitely recognizable in the phase-folded plot of the YSD case, with a much larger transit depth than the previous case and indeed it is correctly detected by TLS. This suggests that the algorithm efficiency could vary according to the number of quarters in a way that is not straightforward, thus it is recommended that two algorithms be used, instead of just one, when the dataset expands beyond one single quarter.

Finally, considering the entire LCs and the combined results of the three algorithms, we find that merging together the results of YSD and Huber spline increases the final detection efficiency close to 100% for all planetary types, and specifically in Table 1 the final recovery percentages are reported.

We also find that adding the Earth-size planets recovered by N&L to those detected by YSD and Huber spline only slightly improves the final percentage of recovered transits, up to 99.0%. Given that N&L is much more computationally expensive than the other two algorithms, we do not recommend it to filter PLATO LCs encompassing several quarters.

Fig. 13 Injection-recovery tests on the sample of active stars. Injected scaled planetary radius as a function of the injected orbital period: blue dots are the recovered planets, whereas red stars are the missed planets after filtering the LCs with N&L (top) and Huber spline (bottom).

Fig. 14 Results of the planetary injection-recovery tests after filtering the LCs of the quiet sample with N&L. Upper panel: recovered orbital period (Porb) as a function of the corresponding injected Porb. Lower panel: recovered scaled planetary radius (Rp/R*) as a function of the corresponding injected Rp/R*, for the three planetary signals, colored as in the legend. Black solid line marks points where all the data should lay if the recovered values were exactly the same as the injected ones.

Fig. 15 Raw flux (upper panel), detrended flux (middle panel), and phase-folded flux (lower panel) of the LC #58 of the sample of active stars with injected a transit signal with Porb = 32.83 days and Rp/R* = 0.015. Data points are shown in black, whereas different colors represent the photometric variability model inferred from algorithm: N&L (orange) and YSD (red). The transit is undetected in the former case, while it is detected with SDE = 25 in the latter case. Blue vertical lines mark the mid-transit time.

5.4 Peculiar cases

We investigated two peculiar cases that represent a small fraction of the existing planets but that are important from the planetary evolution point of view. The first one are eccentric hot Jupiters (described in Sect. 5.4.1) and the second one is represented by those planets with an orbital period equal to the rotational period of their host star (reported in Sect. 5.4.2). We performed this analysis only on the LCs of the active sample and considering individual quarters, as in Sect. 5.1.

5.4.1 Eccentric hot Jupiters

Hot Jupiters are large planets on extremely close-in orbits with periods shorter than 10 days. The formation pathways for these highly irradiated planets are still under investigation and among them, there are mainly two migration mechanisms: the so-called smooth migration and violent migration (Bailey & Batygin 2018). The latter scenario involves a series of high eccentricity trajectories as a consequence of gravitational interactions with the other bodies of the system, such as planet-planet scattering (Beaugé & Nesvorný 2012) or the Kozai–Lidov mechanism (Shevchenko 2017), until the planet is tidally captured onto a close-in circular orbit around the host star. Studying the detection efficiency of the filtering algorithms for these kind of transit signals is particularly interesting because determining if and how the frequency and the properties (such as the orbital eccentricity) of hot Jupiters vary as a function of age, allows us to understand which of the proposed evolutionary mechanisms actually occur. Therefore, we synthesized a set of 800 transit signals with the following properties: R_p/R_* = 0.1, argument of pericenter (ω) from a uniform distribution in the interval [0–360]°, moving on orbits inclined to 90 degrees, with an eccentricity of 0.2, 0.4, 0.6, or 0.8 and an orbital period in the range of 1–10 days (again sampled from a uniform distribution). Indeed, eccentricity changes the transit shape, depending in turn on the argument of pericenter, as shown in Fig. 16, and we wanted to verify that also this kind of exoplanets could be detected.

We find that eccentricity does not really influence the detection efficiency; in fact, 100% of injected signals were recovered after the filtering of all the selected algorithms. This result was partially expected since simulated transits belonging to this set represented all large planets, which are generally more easily detectable. Besides the orbital periods and mid-transit times, which are also the reference parameters for the recovery of injected signals, the other recovered orbital and planetary parameters are often not in agreement with the corresponding injected values. This happens especially for planets in very eccentric orbits, when the transit fit performed with the TLS default template (with ecc = 0) does not perfectly match the data-points. However, we verified that the correct orbital and planetary parameters can be retrieved by performing a MCMC fit. Moreover, in the VARLET case the algorithm reduces the transit depth of large planets, as already noticed in the previous tests. As a consequence, the SDE is decreased down to a maximum value of about 53, while the maximum SDE value of the recovered planets is around 74 for the other algorithms.

5.4.2 Stellar rotation-planetary orbit commensurability

The second peculiar case is characterized by those planets which orbit with a period exactly equal to their host star rotation period, namely P_orb,planet = P_rot,star. Several stars are observed hosting close exoplanets (mostly giant planets) rotating in an orbit with a period equal to an integer multiple of the stellar rotation period (Walker et al. 2008; Pagano et al. 2009; Szabó et al. 2012; Béky et al. 2014). This is reflected in the photometric observations of the host star which exhibits variations synchronous to the orbital period of the corresponding planet. Studying exoplanetary systems with a stellar rotation-planetary orbit commensurability is important in order to better understand how interactions between the host star and the planet work. Nevertheless, it is challenging to find transiting exoplanets whose orbital period closely matches the rotation period of the star or is commensurable to it, since the transit signal feature may be indistinguishable from the rotational variability of the host star.

This case was tested in this work only on the first quarter of each LC of the active sample. We generated a library of 300 different template transits assuming a circular orbit with an impact parameter b uniformly sampled in [0,1], and with a planet-to-star radius randomly drawn in the range R_p/R_* ∈ [0.01–0.1], dividing the set into three planetary types depending on R_p/R_*, as previously done (see Sect. 5): Jupiter-, Neptune-, and Earth-size planets. Specifically, we produced 100 transit templates with planet-to-star radius in [0.01––0.03] (i.e., Earth-size), 100 with R_p/R_* ∈ [0.03–0.05] (i.e., Neptune-size), and 100 with R_p/R_* ∈ [0.05–0.1] (i.e., Jupiter-size). For all these planets, we chose the orbital period to be exactly equal to the rotation period of their host star.

The percentage of recovered transits of different planetary type is shown in Fig. 17. Since the simulated LCs of the active sample have P_rot ≤ 10 days, there are many injected transits representing ultra-short-period planets, meaning a category of exoplanets defined by having orbital period shorter than a single day. This kind of planet is not so rare, they are about as common as hot Jupiters. Indeed it is estimated that about one out of 200 Sun-like stars has an ultra-short-period orbiting planet (Winn et al. 2018).

From Fig. 17, it can be noticed that, generally, all the filtering algorithms (except N&L) show a good performance for Jupiter-and Neptune-size planets, with recovery efficiency higher than 98%, while in the case of Earth-size planet only biweight and VARLET reach this threshold. Indeed, N&L under-performs with respect to the other algorithms, recovering only 89.33% of the total planets. The missed planets are those with an orbital period shorter than 2 days, namely, those for which also the LOCoR filtering method is used (as shown in Fig. 18), where the scaled planetary radius as a function of the injected orbital period is shown. This result was partially expected since also Rizzuto et al. (2017, authors of this code) noticed that the LOCoR algorithm, which is applied to the most rapidly rotating stars (i.e., P_rot ≤ 2 days), tends to fail in modeling the photometric variability without affecting the in-transit points when the planetary orbital period closely matches the stellar rotation period.

The general-purpose algorithms perform very well for this kind of planets, specifically, the biweight with a window-size of 1 day shows the highest detection efficiency, recovering up to 99% of injected transits of small planets. The Huber spline shows a very high recovery efficiency as well, although it is lower than that of VARLET.

Among the custom-built algorithms, VARLET shows the highest recovery percentage, detecting 99.33% of the total injected transits, and, in particular, 98% of the small planets, a surprising result when compared to those detailed in Sect. 5.1. Although this algorithm achieves a detection efficiency of 100% for larger planets, it strongly reduces their transit depth, as already observed, recovering them with a low SDE of about ~54. For larger planets similar performances are achieved by YSD, where planets are recovered with a higher SDE (up to 90) than VARLET, although some are missed.

In Fig. 19, the SDE boxplot summarizing the properties of the SDE distributions is shown. With regards to the SDE values of YSD and biweight (yellow and orange bars, respectively) they show very similar results, even if the biweight median SDE is slightly larger than the median of YSD. Nevertheless, the highest median SDE is obtained by the Huber spline algorithm, which, however exhibits slightly worse performances than the other two algorithms.

Fig. 16 Illustration of transit shape variation at different orbital eccentricities and for two different argument of pericenter: ω = 90° (top) and ω = 230° (bottom).

Fig. 17 Recovery fraction of injected transit signals of planets with orbital period equal to the stellar rotation period in simulated PLATO stellar LCs of a sample of active stars.

Fig. 18 Injection-recovery tests on the sample of active stars for planetary transit signals with Porb,planet = Prot,star. Injected scaled planetary radius as a function of the injected orbital period: blue dots are the recovered planets, whereas red stars are the missed planets after filtering the LCs with Notch & LOCoR.

Fig. 19 Boxplot of the SDE for the algorithms listed on the y-axis, tested on the active sample with injected transit signals with Porb,planet = Prot,star. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines show the median. The black circles in the VARLET case represent data that extend beyond the whiskers.

6 Discussions and conclusions

The idea behind this work was inspired by the PLATO Working Package (WP) 111 000 (Coordination of Tools for Lightcurve Filtering), a sub-group of the WP 11 that deals with PLATO’s Exoplanet Science. In preparation for this space mission, it was necessary to determine which is the best-performing algorithm for the filtering of stellar activity in PLATO LCs, for the subsequent detection of planetary transit signals. The results presented in this work suggest that the best strategy to adopt when analysing PLATO LCs should be as follows:

• Using the N&L algorithm to filter activity on LCs obtained every three months (i.e., each time that a new quarter of data is added to the previous time series), since it turns out to be the best-performing algorithm on partial time series for both the quiet and active sample. Alternatively, similar results can be achieved by combining the Huber spline and the biweight algorithms.

• After two years, at the end of the observations, that is, once that the LCs of the targets are completed, using both Huber spline and YSD to filter stellar activity on the entire LC (i.e., all the 8 quarters). Indeed, from this analysis, it appears that the combination of these two algorithms succeeds in detecting small planets that are missed by N&L on long temporal baselines, thus achieving the highest possible efficiency.

• For stars with a rotation period shorter than 2 days²⁴, it is recommended to use the biweight algorithm on the entire LCs, since it shows the highest recovery efficiency and the highest SDE values in detecting small exoplanets with an orbital period that closely matches the stellar rotation period. This kind of transit signals are often missed by N&L. We note that in this case, VARLET is the second best-performing algorithm despite the lower average SDE, in opposition to the general trend observed in the other cases.

Among the tested algorithms, the GPs showed poor performances compared to the others. However, we would like to emphasize that we investigated Gaussian Processes using only a specific kernel (the Celerite rotation term implemented in celerite2) and that the results could significantly change with a different kernel, such as a Matern 3/2 term or a squared exponential. It is well known that Gaussian processes tend to absorb the transit signal, by modeling it as part of the photometric variability, and that this problem can be alleviated by simultaneously fitting for the planetary transit. We also note that the biweight algorithm can achieve better performances by fine-tuning the window width to three times the transit duration, as found by H19, although in our case, there is no prior knowledge of the orbital parameters of the planet. In our analysis, we investigated the biweight method using a fixed window of 1.0 day as the best approximation given the expected maximum transit duration of our injected transit models; however, possible performance gains could be achieved with a different value or by applying an adaptive window width.

In a recent paper, Nardiello et al. (2021) investigated the stellar age versus planetary radius distribution by using the results of different works on both open clusters and young stellar associations, with well-constrained stellar ages, containing stars hosting both confirmed and candidate exoplanets. This distribution is shown in Fig. 7 of Nardiello et al. (2021), where it is clearly visible that small exoplanets, namely, objects with R_p ≤ 4R_Earth, are all concentrated at ages older than 100–200 Myr, whereas larger planets exhibit a broader distribution. This concentration of small objects around older stars might be due to either planetary formation and evolution processes or to an observational bias, since it has been particularly difficult to detect transits of Earth and Super-Earth-size planets in the LC of young active stars to date. Until a systematic analysis over all ranges of planetary radius is performed, it is difficult to draw any theoretical conclusion with certainty. To date, such an analysis cannot be executed with current astrophysical instrumentation, but it may be possible thanks to PLATO in the coming years. Indeed, according to the results of this work, it seems that when turning from TESS to PLATO LCs of young active stars, the efficiency of the filtering algorithms for the planetary transit detection increases from about 40% up to nearly 99%. The injection-retrieval experiments performed on PLATO simulated LCs have shown that PLATO will be able to detect a large variety of exoplanets with many different characteristics, in terms of size, orbital periods, and eccentricity. In particular, they have demonstrated that detecting Earth-size exoplanets with PLATO will be possible with both general-purpose and custom-built filtering algorithms, even in stellar LCs of young active stars. Therefore, in the forthcoming years only PLATO can tell us if the lack of small exoplanets around young stars is really related to an observational bias.

Even if the injection-recovery tests performed in this work included several different kinds of planets, we have not explored yet the entire parameter space of the observed exoplanets. Considering the distribution of planets detected by the transit technique, we are far from completion. Other peculiar cases that could be tested and simulated in future works include, for example, transits showing timing variations (TTV; Holman & Murray 2005; Nascimbeni et al. 2011; Malavolta et al. 2017), monotransits (Cooke et al. 2020), multi-planetary systems, or small planets with orbital periods longer than 40 days, that is, temperate terrestrial planets such as those around the dwarf star TRAPPIST-1 (Gillon et al. 2017). Indeed, PLATO’s top-level science requirements are focused on habitable planets, therefore it would be interesting to test the algorithms also on more temperate Earth-like planets. Efforts to compute the expected yield of Earth analogs by PLATO are already ongoing, for example, by Heller et al. (2022).

Furthermore, in future works, other filtering algorithms can be tested in order to analyze their efficiency on PLATO LCs and to compare the results with those presented in this paper. For example, the Kepler Science Data Processing pipeline²⁵ (Jenkins et al. 2010) that is now publicly available, has not been tested on active stars yet, but it has demonstrated very good results on quiet stars, so it would be interesting to check its performance on the active sample of PLATO LCs.

Finally, we note that for our experiment, we employed the TLS algorithm only, assuming that a transit-like shape could be more efficient than a box-like shape in detecting the planet. However, the deformation of the transit feature introduced by a filtering algorithm may invalidate this assumption (see Fig. 10). Future works related to the planetary yield of PLATO should consider multiple filtering algorithms as well as multiple detection methods. It is also worth noting that the full PLATO pipeline is now under development, as well as extended simulations including the effects of contaminants in the LC, while a complete and robust assessment of the in-flight performance will be available only after the launch. Although we could have removed the long-term instrumental systematics; for instance, with the use of a Gaussian Process, we preferred to use the LC as they are produced by the PLATO Solar-like Light-curve Simulator and test the worst-case scenario in which an intermediate data product has to be analyzed. As our analysis is purely differential, that is, we are not trying to determine PLATO’s yield for a specific class of planets, our results should hold as long as the comparison is performed in a self-consistent way and the main properties of the simulations resemble those of the final LCs. When simulated final PLATO data products become available, we plan to perform new injection-retrieval tests, including both peculiar cases not considered in this work, such as mono-transits or long-period planets, as well as custom-built filtering algorithms that will be presented to the community.

Acknowledgements

This work presents results from the European Space Agency (ESA) space mission PLATO. The PLATO payload, the PLATO Ground Segment and PLATO data processing are joint developments of ESA and the PLATO Mission Consortium (PMC). Funding for the PMC is provided at national levels, in particular by countries participating in the PLATO Multilateral Agreement (Austria, Belgium, Czech Republic, Denmark, France, Germany, Italy, The Netherlands, Portugal, Spain, Sweden, Switzerland, Norway, and UK) and institutions from Brazil. Members of the PLATO Consortium can be found at https://platomission.com/. The ESA PLATO mission website is https://www.cosmos.esa.int/plato. We thank the teams working for PLATO for all their work. L.M., I.P., G.P., and S.D. acknowledge support from PLATO ASI-INAF agreements no. 2015-019-R0-2015 and no. 2015-019-R.1-2018.

Appendix A Orbital and planetary parameters from TLS and MCMC fit of LC #17

An example of a grazing transit when the TLS default template does not correctly fit the data-points (and, consequently, it fails to retrieve most of the orbital and planetary parameters) is given in Fig. A.2. This transit represents a large planet with a R_p/R_* ~ 0.09, whereas the TLS recovered a planetary signal corresponding to a much smaller body of only R_p/R_* ~ 0.037, as shown in Fig. A.2 (left), where the phase-folded data and the residuals of the transit fit performed with TLS are shown in the upper and lower panel, respectively. In Fig. A.1 the filtering performed by YSD on the raw flux is shown. Even when performing a fit with the TLS grazing template, the inferred scaled planetary radius is 0.043, slightly improving with respect to the TLS default template fit but still underestimating the injected value. Therefore, a Markov chain Monte Carlo (MCMC) analysis has been run in order to fit the correct transit model to the light curve and thus achieve a more accurate planet characterization. For this purpose, the PyORBIT²⁶ code (Malavolta 2016) was used to perform a MCMC analysis with 16 walkers and 500000 steps (the burn-in is after 20000 steps), initialized with a normal distribution around stellar radius and density. Moreover, uniform priors for the other parameters have been applied, such as the period and the time of transit, centered on the values obtained with TLS, within their boundaries. We further assumed, as priors on the limb darkening coefficients u₁ and u₂, a Gaussian distribution centered on the values given by tables describing the simulated LCs, provided by LSWG. The relevant physical parameters obtained from the posterior samples of the best-fit results are reported in Table A.1.

Uncertainties on the orbital and planetary parameters were inferred from their a posteriori distributions obtained from the MCMC analysis. The transit fit and its residuals are shown in Fig. A.2 (right). It is possible to notice that the transit fit significantly improves by performing the MCMC analysis, recovering parameter values consistent with the injected ones. This suggests that the scaled planetary radius obtained with TLS must be considered carefully, since the real transit shape could significantly differ from the TLS template, thus resulting in a systematic error which gives an under-estimate of R_p/R_* – consequently affecting also all the other parameters related to R_p/R_*. On the other hand, the estimates of the orbital period and the time of transit recovered from the TLS transit search result to be very precise and generally reliable.

Finally, we emphasize that the TLS algorithm has not been designed for transit characterization, but only for transit detection. Indeed, an MCMC analysis generally needs to be performed in order to characterize real planetary transits.

Appendix B Stellar parameters of the simulated LCs

The simulated light curves of both the quiet and active sample are provided as supplementary material. In Fig. B.1, Fig. B.2, Fig. B.3, and Fig. B.4, a selection of four simulated LCs from the quiet sample, before and after filtering is shown. A large planet (i.e., R_p/R_* ≥ 0.05) is injected, as explained in Sect. 5. In these figures, six panels are displayed, representing the filtering performed by the selected algorithms. Specifically, in the left panels, we show the raw flux with the superimposed model, whereas in the right panels, we show the resulting filtered flux. Different colors are adopted for the model produced by different algorithms. In particular: Yellow for N&L, red for YSD, fuchsia for GPs, violet for VARLET, blue for biweight, and green for the Huber spline. As a comparison, also a selection of 4 LCs from the active sample is shown in Fig. B.5, B.6, B.7, and B.8. Moreover, the full set of input parameters for the 100 simulated PLATO LCs of the quiet and active sample is reported in Table B.1 and Table B.2, respectively.

Fig. A.1 Filtering by YSD of the first quarter of simulated LC #17 of the quiet sample with injected a large planet. In blue, we show the raw flux, in black the filtering model, and in red the corrected flux are shown, as in the legend.

Fig. A.2 Transit fit of the first quarter of LC #17 with a large injected planet. Top panel: Best fit of the transit light curve is represented by the solid black line, whereas the blue points are the PLATO normalized data, as a function of the time from mid-transit. Bottom panel: Residuals (O-C: observed data-calculated data) are shown as a function of the time from mid-transit. Left: Transit model computed with TLS default template. Right: Transit model computed after an MCMC analysis.

Fig. B.1 Filtering of quarter 7 of LC #2 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.2 Filtering of quarter 6 of LC #24 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.3 Filtering of quarter 5 of LC #30 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.4 Filtering of quarter 4 of LC #47 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.5 Filtering of quarter 1 of LC #17 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.6 Filtering of quarter 5 of LC #29 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.7 Filtering of quarter 2 of LC #31 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Fig. B.8 Filtering of quarter 1 of LC #97 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.

Appendix C Comparison of running times

In Table C.1, we report the running times of the filtering algorithms analyzed in this work, we computed over 20 PLATO simulated LCs of the active sample encompassing one single quarter (including 12486 data points). The analysis was carried out on an outdated desktop machine hosting two Intel(R) Xeon(R) E5-2620 at 2.00GHz CPUs for a total of 24 threads and a 2TB hard drive spinning at 5400 RPM. We did not make any attempt to optimize the (extremely heterogeneous) preexisting implementations of the algorithms or to allow for multiprocessing analysis if not already implemented, as this task was beyond the scope of this work. For these reasons, we did not consider the execution time a relevant factor in our comparisons.

References

All Tables

All Figures

	Fig. 1 PLATO simulated LC #0 of the sample of quiet stars (top) and active stars (bottom). Different colors represent different quarters. The systematics in the simulated LCs represent a worst-case scenario, as their correction will be addressed in the final PLATO data products.
In the text

Fig. 2 Recovery fraction of injected transit signals similar to those of Hippke et al. (2019) injected into simulated PLATO stellar LCs of quiet (second column) and active (third column) stars, compared with the recovery fraction of similar signals injected into a sample of active stars observed by TESS (first column) and analyzed by H19. The biweight with a window of 1.0 day was not tested by H19 on TESS active LCs.

In the text

Fig. 3 Raw flux (left panel) and detrended flux (right panel) of the 3rd quarter of LC #9 of the sample of active stars with an injected transit signal with Porb = 12.60 days and Rp/R*=0.056. Data-points are shown in black, whereas different colors represent the photometric variability model inferred from algorithm: Median (red), biweight with a window-size of 0.25 days (dark blue) and 1.0 day (light blue), Huber spline (green) and Lowess (orange).

In the text

Fig. 4 Boxplot of the signal detection efficiency for the general-purpose algorithms from Wōtan tested on the active sample with injected transit signals with Rp ~ 0.5 RJ and Porb ∈ [1–15] days, moving on circular orbits. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines the median. The black circles in the biweight (1.0 day) and the lowess (1.0 day) cases represent data that extend beyond the whiskers.

In the text

	Fig. 5 First quarter of the PLATO simulated LC #0 of the quiet (left) and active (right) sample, with an injected Jupiter-size planet. In black, we show the normalized flux without transits and in red, the transit model generated with batman.
In the text

	Fig. 6 First quarter of the PLATO simulated LC #0 of the quiet (left) and active (right) sample, with an injected Neptune-size planet. In black, we show the normalized flux without transits and in blue, the transit model generated with batman.
In the text

	Fig. 7 First quarter of the PLATO simulated LC #0 of the quiet (left) and active (right) sample, with an injected Earth-size planet. In black, we show the normalized flux without transits and in green, the transit model generated with batman. In the active case, the transits are hard to distinguish, therefore a zoom-in on the first transit is shown. The second transit occurs at about 70 days.
In the text

	Fig. 8 Recovery fraction of injected transit signals in simulated PLATO stellar LCs of the sample of quiet stars.
In the text

	Fig. 9 Recovery fraction of injected transit signals in simulated PLATO stellar LCs of a sample of active stars.
In the text

	Fig. 10 Raw flux (top) and detrended flux (bottom) of the first quarter of LC #22 of the quiet sample, zoomed in on the transit points of the injected Earth-size planet. Data points are shown in black, whereas different colors represent the photometric variability model inferred from different algorithms: N&L (orange), YSD (red), GPs (fuchsia), and VARLET (violet).
In the text

	Fig. 11 Boxplot of the SDE for the algorithms listed on the y-axis, tested on the quiet sample with injected transit signals as described in Sect. 5.1. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines the median. The black circles in the VARLET case represent data that extend beyond the whiskers.
In the text

	Fig. 12 Boxplot of the SDE for the algorithms listed on the y-axis, tested on the active sample with injected transit signals as described in Sect. 5.1. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines the median. The black circles in the VARLET case represent data that extend beyond the whiskers.
In the text

	Fig. 13 Injection-recovery tests on the sample of active stars. Injected scaled planetary radius as a function of the injected orbital period: blue dots are the recovered planets, whereas red stars are the missed planets after filtering the LCs with N&L (top) and Huber spline (bottom).
In the text

Fig. 14 Results of the planetary injection-recovery tests after filtering the LCs of the quiet sample with N&L. Upper panel: recovered orbital period (Porb) as a function of the corresponding injected Porb. Lower panel: recovered scaled planetary radius (Rp/R*) as a function of the corresponding injected Rp/R*, for the three planetary signals, colored as in the legend. Black solid line marks points where all the data should lay if the recovered values were exactly the same as the injected ones.

In the text

Fig. 15 Raw flux (upper panel), detrended flux (middle panel), and phase-folded flux (lower panel) of the LC #58 of the sample of active stars with injected a transit signal with Porb = 32.83 days and Rp/R* = 0.015. Data points are shown in black, whereas different colors represent the photometric variability model inferred from algorithm: N&L (orange) and YSD (red). The transit is undetected in the former case, while it is detected with SDE = 25 in the latter case. Blue vertical lines mark the mid-transit time.

In the text

	Fig. 16 Illustration of transit shape variation at different orbital eccentricities and for two different argument of pericenter: ω = 90° (top) and ω = 230° (bottom).
In the text

	Fig. 17 Recovery fraction of injected transit signals of planets with orbital period equal to the stellar rotation period in simulated PLATO stellar LCs of a sample of active stars.
In the text

	Fig. 18 Injection-recovery tests on the sample of active stars for planetary transit signals with Porb,planet = Prot,star. Injected scaled planetary radius as a function of the injected orbital period: blue dots are the recovered planets, whereas red stars are the missed planets after filtering the LCs with Notch & LOCoR.
In the text

Fig. 19 Boxplot of the SDE for the algorithms listed on the y-axis, tested on the active sample with injected transit signals with Porb,planet = Prot,star. Boxes cover the lower to upper quartiles; whiskers show the 10 and 90 percentiles. Dashed black lines indicate the mean SDE, whereas red lines show the median. The black circles in the VARLET case represent data that extend beyond the whiskers.

In the text

	Fig. A.1 Filtering by YSD of the first quarter of simulated LC #17 of the quiet sample with injected a large planet. In blue, we show the raw flux, in black the filtering model, and in red the corrected flux are shown, as in the legend.
In the text

Fig. A.2 Transit fit of the first quarter of LC #17 with a large injected planet. Top panel: Best fit of the transit light curve is represented by the solid black line, whereas the blue points are the PLATO normalized data, as a function of the time from mid-transit. Bottom panel: Residuals (O-C: observed data-calculated data) are shown as a function of the time from mid-transit. Left: Transit model computed with TLS default template. Right: Transit model computed after an MCMC analysis.

In the text

	Fig. B.1 Filtering of quarter 7 of LC #2 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.2 Filtering of quarter 6 of LC #24 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.3 Filtering of quarter 5 of LC #30 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.4 Filtering of quarter 4 of LC #47 of the quiet sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.5 Filtering of quarter 1 of LC #17 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.6 Filtering of quarter 5 of LC #29 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.7 Filtering of quarter 2 of LC #31 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text

	Fig. B.8 Filtering of quarter 1 of LC #97 of the active sample with an injected large planet. The filtering models are overplotted in different colors and the resulting detrended flux is shown on the right panels.
In the text