Accounting for stellar activity signals in radial-velocity data by using Change Point Detection techniques

Active regions on the photosphere of a star have been the major obstacle for detecting Earth-like exoplanets using the radial velocity (RV) method. A commonly employed solution for addressing stellar activity is to assume a linear relationship between the RV observations and the activity indicators along the entire time series, and then remove the estimated contribution of activity from the variation in RV data (overall correction method). However, since active regions evolve on the photosphere over time, correlations between the RV observations and the activity indicators will correspondingly be anisotropic. We present an approach that recognizes the RV locations where the correlations between the RV and the activity indicators significantly change in order to better account for variations in RV caused by stellar activity. The proposed approach uses a general family of statistical breakpoint methods, often referred to as change point detection (CPD) algorithms; several implementations of which are available in R and python. A thorough comparison is made between the breakpoint-based approach and the overall correction method. To ensure wide representativity, we use measurements from real stars that have different levels of stellar activity and whose spectra have different signal-to-noise ratios. When the corrections for stellar activity are applied separately to each temporal segment identified by the breakpoint method, the corresponding residuals in the RV time series are typically much smaller than those obtained by the overall correction method. Consequently, the generalized Lomb-Scargle periodogram contains a smaller number of peaks caused by active regions. The CPD algorithm is particularly effective when focusing on active stars with long time series, such as alpha Cen B.


Introduction
The radial velocity (RV) method (e.g., Mayor & Queloz 1995;Lovis & Fischer 2010;Hatzes 2016) is one of the most successful techniques for detecting extrasolar planets (e.g., Fischer et al. 2016). Nevertheless, perturbations in the RV data caused by different types of stellar signals such as stellar oscillations, granulations, and the presence of active photospheric regions continuously plague the exoplanet community (e.g., Saar & Donahue 1997;Queloz et al. 2001;Lindegren & Dravins 2003;Desort et al. 2007;Meunier et al. 2010Meunier et al. , 2017Dumusque et al. 2011a;Dumusque 2016Dumusque , 2018. In particular, stellar activity in the form of spots and faculae represents the main limitation in the full characterization of Earth-like exoplanets in RV surveys (e.g., Saar & Donahue 1997;Meunier et al. 2010;Dumusque et al. 2014;Dumusque 2018). New generation spectrographs such as the EXtreme PREcision Spectrometer (EXPRES, Jurgenson et al. 2016) and the Echelle SPectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO, Pepe et al. 2014) have recently been constructed to solve the issue of instrumental signal, and therefore to provide data only affected by stellar and planetary signals. However, disentangling the perturbations due to stellar activity from the RV variations caused by small size exoplanets remains the most crucial challenge (e.g., Davis et al. 2017;Dumusque 2018;Reiners et al. 2018;Simola et al. 2019), as the RV variations induced by active regions are an order of magnitude larger than those expected from Earth-like exoplanets (e.g., Dumusque 2018).
Efforts have been made to model the signals caused by different sources of stellar variability within the RV time series (e.g., A127, page 1 of 29 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication. Tuomi et al. 2013;Rajpaul et al. 2015;Davis et al. 2017;Simola et al. 2019). Several solutions have been successfully proposed in order to deal with stellar oscillations and granulation phenomena, such as: calculating stellar evolution sequences (e.g., Christensen-Dalsgaard et al. 1995); fitting a two-level structure tracking (TST) algorithm based on a two-level representation of granulation (Del Moro 2004); using daytime spectra of the Sun in order to measure the solar oscillations (e.g., Kjeldsen et al. 2008;Lefebvre et al. 2008); and characterizing the statistical properties of magnetic activity cycles focusing on HARPS observations (e.g., Pepe et al. 2011;Dumusque et al. 2011b). However, properly modeling the other sources of stellar activity remains extremely challenging (e.g., Nava et al. 2019). In the present work, we deal with the cross correlation function (CCF) that is derived from the stellar spectrum (e.g., Hatzes 1996;Hatzes & Cochran 2000;Fiorenzano et al. 2005). As it is well known, the CCF barycenter estimates the RV of the star. The asymmetry and the full width at half maximum (FWHM) of the CCF give a strong indication for stellar activity, meaning that variations in RV are caused by active regions rather than by an exoplanet (e.g., Hatzes 1996;Queloz et al. 2001;Boisse et al. 2011;Figueira et al. 2013;Simola et al. 2019). Several solutions have been successfully proposed for mitigating stellar activity perturbations when working with RV measurements, including: decorrelating RV data against activity indicators such as log R HK (e.g., Wilson 1968;Noyes et al. 1984) or Hα (e.g., Robertson et al. 2014); modeling stellar activity by fitting Gaussian processes (GPs, Rasmussen & Williams 2005;Haywood et al. 2014;Rajpaul et al. 2015); or moving averages (e.g., Tuomi et al. 2013) to the RV data. A common statistic employed for identifying changes in the shape of the CCF is the bisector span (e.g., Hatzes 1996;Queloz et al. 2001).
As already pointed out, disentangling the signal caused by stellar active regions and the signal due to an Earth-like planet is extremely challenging, but some differences in the characteristics of the two signals may be of help. Active regions produced by spots and faculae evolve on the star's photosphere and can generate RV variations that spread from a few days to several weeks, depending on the rotation period of the star and the intensity of the magnetic cycle (e.g., Noyes et al. 1984;Dumusque et al. 2014;Davis et al. 2017;Dumusque 2018;Nava et al. 2019). While exoplanets produce a Doppler-shift on the CCF without modifying its shape or its width, active regions produce variations in both the asymmetry and the FWHM of the CCF, whose effect is a variation in the barycenter of the CCF, and therefore in the estimated median of the RV CCF (see e.g., Hatzes 1996;Queloz et al. 2001;Boisse et al. 2011;Figueira et al. 2013;Dumusque et al. 2017;Simola et al. 2019; Thompson et al. 2020), hereafter indicated as RV. We refer the reader to Table 1 for the main RV-related symbols we adopt throughout paper. Moreover, while an exoplanet produces a persistent Keplerian signal, the signal produced by active regions is not persistent as it waxes, wanes, and changes as a function of time (e.g., Fischer et al. 2016;Davis et al. 2017;Dumusque 2018;Thompson et al. 2020).
By using the stellar activity parameters estimated from the study of the CCF, a linear model is often proposed in order to decorrelate the RV data against the RV variations due to stellar activity (e.g., Dumusque et al. 2017;Feng et al. 2017a;Simola et al. 2019). Rather than using a normal fit to the CCF and then calculating the bisector span, Simola et al. (2019) used a skew normal (SN) fit to the CCF. The SN distribution allows a parameter to be specified, hereafter γ, which describes the asymmetry of the CCF (e.g., Simola et al. 2019;Adcock & Azzalini 2020). Table 1. Legend of the RV-related symbols adopted in this paper.

RV Radial velocity RV
Median of the RV CCF based on an SN fit RV * activity β 0 + β 1 A + β 2 γ + β 3 FWHM SN  Notes. See the text for further details. Simola et al. (2019) suggest using the median RV as the barycenter of the SN fit to the CCF. Given the benefits of fitting an SN to the CCF as emphasized by Simola et al. (2019), the model we employed in the analyses of this work is defined as follows: where β 0 is the coefficient corresponding to the intercept, A is the contrast parameter of the CCF (Fig. 2 in Dumusque et al. 2014), γ is the aforementioned CCF asymmetry parameter, FWHM SN is the FWHM of the CCF fitted with the SN, and is the random error having a multivariate normal distribution with vector of means equal to 0 and variance-covariance matrix equal to σ 2 I (with I being the identity matrix).
The temporal evolution of active regions is reflected by the change of the activity indicators obtained from the CCF (i.e., the already defined A, γ, and FWHM SN ) and depends on the rotation period of the star and on its magnetic cycle (e.g., DeWarf et al. 2010;Dumusque et al. 2011b;Borgniet et al. 2015). It seems therefore reasonable to assume that stellar activity is not stationary, but rather "piecewise stationary". By the term "piecewise stationary", we mean that the stellar activity does not change significantly within certain properly selected temporal segments. If we were able to detect the bounds of those segments, we could properly split the RV time series and correct for stellar activity by applying Eq. (1) to each segment. This represents the core of the change point detection (CPD) methods to be compared with the overall correction (oc) method, where the correction model based on Eq. (1) is applied over the entire RV time series.
In this paper, we propose the breakpoint (bp) method using a CPD technique to correct for variations within an RV time series that are induced by active stellar regions. The bp method is compared to the oc method using four different stars, some of which have known exoplanets. We demonstrate that the bp method is better able to correct the RV time series variations, suggesting that the class of CPD methods may be helpful for detecting low-RV planetary signals such as the ones induced by Earth-like exoplanets.
The paper is organized as follows. In Sect. 2, we introduce the CPD methods, which constitute the statistical framework we used to perform our analyses. In Sect. 3, we compare the performances of the CPD method in use and the oc method, both using real observational data and developing a proper simulation study to quantify the threshold of detection of exoplanets. The discussion of the results and our conclusions are outlined in Sects. 4 and 5, respectively. A127, page 2 of 29 U. Simola et al.: CPD exoplanets

Change point detection methods
CPD methods are widely used when the goal is to find changes and variations in a time series (e.g., Truong et al. 2020). The presence of change points in a time series is a strong indication that the data generating process has changed (e.g., van den Burg & Williams 2020). The first CPD method (the so-called "intercept-only" method) was originally introduced by Page (1954) in order to identify changes in the mean of an industrial quality control variable. CPD methods have since become a reliable and widely used solution in bioinformatics (e.g., Guédon 2013;Hocking et al. 2013;Truong et al. 2018), climatology (e.g., Reeves et al. 2007;Verbesselt et al. 2010;Maidstone 2016), financial analyses (e.g., Lavielle et al. 2006;Frick et al. 2014), medicine (e.g., Liu et al. 2018), network data traffic analysis (e.g., Lévy-Leduc 2009;Lung-Yut-Fong et al. 2011), signal processing analysis (e.g., Lavielle & Teyssiere 2007;Jandhyala et al. 2013;Haynes et al. 2017), and speech processing (e.g., Angelosante & Giannakis 2012).
Following Truong et al. (2020), a CPD method may be classified as either online or offline. The online methods are used when seeking real-time changes in the data generating process (e.g., Adams & MacKay 2007;Sahki et al. 2018). Instead, the so-called offline (referred to also as "a posteriori") methods are based on the CPD algorithms designed to estimate change points in the generative process of a time series when the collecting data process is over (e.g., Truong et al. 2020;van den Burg & Williams 2020). If only a single-parameter change is monitored, then we talk about univariate processes and univariate CPD methods being employed; otherwise we deal with multivariate methods (e.g., Aminikhanghahi & Cook 2017).

The CPD statistical framework
In order to introduce the CPD methods, we start by defining the object of interest as an univariate time series y = y 1,...,T = {y i } T i=1 made of T data points. The y time series is assumed to be "piecewise stationary" (e.g., Chen & Gupta 2011;Truong et al. 2020;van den Burg & Williams 2020), which means that the temporal behavior of y changes at D different and (generally) unknown locations. The goal of the CPD methods is to estimate the set of indices l corresponding to the D locations, l = {l 1 , . . . , l D }, which delimit different stationary regions; by convention l 0 ≡ 1 and l D+1 ≡ T .
From a statistical standpoint, CPD methods fall into the model selection problem, because the overall goal of estimating the vector l of change point locations is equivalent to retrieving the best segmentation of the time series among all the possibilities (e.g., Chen & Gupta 2011;Truong et al. 2020;van den Burg & Williams 2020).
In order to retrieve the best possible segmentation, CPD methods rely on the maximization of the penalized loglikelihood function L P (D; y), which is defined as: where k is the index used to identify the D + 1 segments of y bounded by l 0 , . . . , l D+1 , log L(y [l k−1 :l k −1] ) is the log-likelihood function evaluated at the kth segment, λ is a positive pre-selected constant, and P(T ) is a penalty function that may be added to balance out the goodness-of-fit term given by the log-likelihood function L(D; y) ≡ D+1 k=1 log L(y [l k−1 :l k −1] ) (e.g., Zeileis et al. 2002;Bai & Perron 2003;Truong et al. 2020). After evaluating theD value for which L(D; y) is maximum (i.e., L(D; y) = L(y) = max D L(D; y)), λP(T ) (which is, in general, a function ofL(y), see e.g., Eq. (5)) is subtracted from L(D; y) to obtain L P (D; y). The best partitionl = {l k } for that givenD is the one corresponding to the maximum value of the penalized log-likelihood function, that isL P = max l L P (D, y {l} ). In other words, the different L P (D; y {l} ) values are used to address the model selection problem (i.e., the choice of the best segmentation), andL P identifies the optimal partition of the y time series made ofD + 1 segments. If the total number of change points is assigned a priori and only their locations along y is unknown, no penalty is added in Eq. (2).

The breakpoint method
Given the multiple sources of stellar activity, the most adequate methods to treat the changes in RV time series are the offline multivariate CPD methods, which allow a linear relationship to be expressed between the response variable y and a set of covariates (which may be arranged in a matrix, X), as the one defined in Eq. (1). Using those methods, the set of locations l = l 1 , . . . , l D is estimated not only by evaluating the changes in the mean and in the variance of y, but also by considering the relation between y and X (e.g., Hackl & Westlund 1989;Bai 1997a,b;Bai & Perron 1998Zeileis et al. 2002Zeileis et al. , 2003. In the analyses presented in the following Sections, we use the CPD method proposed by Bai & Perron (2003), known as the breakpoint (bp) method.
In order to introduce the bp method, we start by providing the definition of the multivariate linear regression model made of p covariates for the vector of the observations y, which generalizes the model proposed in Eq. (1). In matrix form, this is: where X is the matrix of covariates, β is the vector of parameters to be determined according to the CPD method in use 1 , and is the already defined random error distributed according to a multivariate normal distribution having vector of means equal to 0 and variance-covariance matrix equal to σ 2 I. Assuming the whole time series y as stationary implies that all the coefficients β = (β 0 , . . . , β p ) are constant over the entire temporal range of observations. Conversely, assuming piecewise stationarity, y can be divided in D + 1 segments over each of which the coefficients are constant, so that there is one set of coefficients per segment, for a total of D + 1 vectors β k (k = 1, . . . , D + 1). With the piecewise stationarity assumption, following Bai & Perron (2003); Zeileis et al. (2003), we may update Eq. (3) as: The index k = 1, . . . , D + 1 indicates the temporal segment over which the vector of the β k coefficients is constant. Bai & Perron (2003) retrieve the simultaneous estimation of multiple breakpoints l 1 , . . . , l D by calculating the residual sum of squares (RSS) of Eq. (4). As is distributed according to a multivariate normal distribution, we recall that the least squares solution proposed by Bai & Perron (2003) is equivalent to the maximum likelihood solution retrieved by maximizing Eq. (2), as long as the same penalty is used in both procedures. In particular, following the implementation of Bai & Perron (2003), λ = 1 and the penalty function P(T ) is given by the Bayesian information criterion (BIC, Schwarz 1978), leading to where the multiplicative term (p + 1)(D + 1) is equal to the overall number of parameters used by the bp method. According to Eq. (1), p + 1 = 4 parameters (i.e., β 0 , . . . , β 3 ) are used inside each segment.
Since D > 1 is unknown, the optimization problem might be computationally challenging. In order to tackle the computational burden of estimatingD, Bai & Perron (2003) used the dynamic programming algorithm originally proposed by Fisher (1958) and extended by Bellman & Roth (1969) and Guthery (1974). The core idea of the dynamic programming approach is to first compute a set of triangular RSS matrices for each D = 1, . . . , D max value (where D max is an integer number tuned by the researcher, indicating the maximum number of breaks allowed in that analysis 2 ). In other words, within a given set of triangular RSS matrices, the number of segments D is fixed, while the starting and ending locations of each segment change. Then, the bp method returns that partition, which minimizes the RSS of Eq. (4) (and hence maximizes L P (D; y) of Eq. (2)) among all the previously created matrices. The main challenge of the dynamic programming approach is given by the computational effort needed to compute the sets of triangular RSS matrices. Further details about the dynamic programming algorithm and the bp method can be found in Bai & Perron (2003); Zeileis et al. (2002Zeileis et al. ( , 2003. In particular, this work uses the version proposed by Zeileis et al. (2003).

Preliminary considerations
To test the performances in modeling and removing stellar activity from RV time series, we carried out several analyses of RV time series to compare the bp method presented in Sect. 2.2 with the oc method that corrects for stellar activity over the entire temporal range considered as a whole segment. We considered real RV observations of four different stars, namely α Cen B, HD 215152, HD 10700, and HD 192310, which were selected because of their high signal-to-noise ratio (S/N) and long-baseline coverage with an intensive observational cadence.
In particular, Rajpaul et al. (2015) analyzed the 459 HARPS RV data points of α Cen B obtained by Dumusque et al. (2012) (a subsample of the RV observations used in this work). After properly modeling the RV time series, they concluded that all the significant variations in the time series are compatible with stellar activity. Therefore, as α Cen B does not likely host any planet, by applying both the bp and oc methods to its RV time series, we may test the performances of the two methods in removing stellar sources of noise and retrieving a "cleaned" time series. Only if a further signal is present will it then be worth wondering whether it is exoplanetary in origin.  Baliunas et al. (1996); (d) Pepe et al. (2011). We considered RV data points derived only from spectra with S /N ≥ 10 at 550 nm. The average S/N values at 550 nm, the number of CCF data points, and the rotation period P rot of the four stars analyzed in this work are listed in Table 2.
All the results were obtained using the statistical software R (R Core Team 2019). In particular, the vector listing the locations of change points within the bp method was estimated by using the R package strucchange 3 . Since the presence of outliers in the data could have led to an exaggerated number of change point locations (e.g., Fearnhead & Rigaill 2019), before applying the bp method and the oc method to the data, we preliminarily rejected all those measurements falling off the range between the 5th and 95th percentile for each considered variable (e.g., Ghosh & Vogt 2012; Pollet & van der Meij 2017).

Analyses of a real star: α Cen B
The RV data of α Cen B consist of T = 17627 CCFs taken from the end of February 2008 to the end of May 2013 by the HARPS spectrograph. Firstly we used the SN fit to the CCF in order to estimate the parameters RV, A, γ, and FWHM SN . Secondly, following Dumusque et al. (2012), all the measurements were preliminarily corrected from the contamination that was induced by the close-in α Cen A. Finally, the outliers were clipped as specified above, which resulted in a final set of T = 16451 measurements.
After cleaning the data, the response variable becomes y = RV, while the matrix of the covariates is X = (1, A, γ, FWHM SN ). Since X also contains the intercept term X 0 ≡ 1, we emphasize again that the parameter vector β is made of p + 1 = 4 elements, where p is the number of covariates.
We applied both the oc and bp methods to the y = RV time series, further assuming as input the linear regression model of Eq. (1) that specifies the relationship between the response variable y and the matrix of covariates X. The bp method foundD = 4 change points (i.e.,D + 1 = 5 piecewise stationary segments), whose locations are summarized in Table 3, where the error bars at the1σ level highlight the precision of the bp method. The variability of the activity indicators inside each segment is extremely small, while the variability between each segment is larger, suggesting indeed that the activity of the star changes where a new segment starts. To highlight the locality, the spread, and the skewness of the data synthesized in Table 3  Notes. For each location, the JD and its corresponding date in the Gregorian Calendar format is displayed. For each segment, the number of data points (# CCFs), the response variable RV, and the estimates of the detrending parameters are shown by reporting their median value and error bars at the 1σ level.   Notes. A 95% confidence interval for each correlation coefficient ρ is specified within parentheses. Notes. We recall that the RV statistic depends on the SN fit, and hence the rms values are the same for both the bp and oc methods, which are applied afterward. The last row shows the BIC coefficients corresponding to the bp and oc models.
difference between contiguous boxplots describing the covariate distributions in contiguous segments visually marks the change in the correlation between RV and the activity indicators. From a quantitative point of view, the Mann-Whitney test (McKnight & Najab 2010) -which assumes that the two samples to be compared are not statistically different as null hypothesis -confirms the differences of the covariates distributions between contiguous segments. In fact, for all the covariates, it gives p-values ∼10 −8 for each pair of contiguous segments. Moreover, we note that A globally decreases, while γ and FWHM SN increase by moving from one segment to the next, which suggests a significant temporal change in the stellar activity level. Overall, Table 3 and Fig. 1 show the characteristics and the quality of the optimal segmentation determined by the bp method for α Cen B. Another instrument used to evaluate the quality of the solution recommended by the bp method is the correlation-pairs plot. Figure 2 displays the correlation plots between RV and the activity indicators when the oc method is applied to the data (top left panel) that is to be compared to the results obtained on each of the five segments found by the bp method (other panels). The oc method is not able to catch all the variations over the entire time series. In fact, the correlation parameters sensibly differ from one segment to the other, suggesting that the piecewise stationary assumption of RV as a function of A, γ, and FWHM SN is reasonable, as it is also highlighted in Table 4.
We can further stress the benefits of using the bp method by comparing ∆RV * bp and ∆RV * oc , which are defined as the residuals obtained using the bp and oc methods, respectively, when the estimated activity level of the star RV * activity is subtracted from the RV time series: The rms of ∆RV * oc vs. ∆RV * bp are 3.48 m s −1 and 1.75 m s −1 respectively, so we reach similar conclusions in preferring the piecewise linear regression strategy, which lowers the rms by 50%. In addition, by comparing the rms of the modeled stellar activity with the rms of the RV signal, it turns out that the bp method explains 89% of the variability of the RV signal in terms of stellar activity, while the oc method is able to model only 40% of the RV signal (see Fig. 3). Since we expect the RV time series to be entirely produced by stellar activity, adopting Eq. (1) in each of the D + 1 segments for the bp method significantly improves the detrending performances.
The detailed rms values of the relevant quantities for both methods are listed in Table 5. In particular, the ∆ BIC between the oc and the bp methods is +22439, thus favoring the bp method according to the BIC minimization criterion. We recall that HARPS is built to obtain RV precision of the order of 1 m s −1 .
A fourth approach used to compare the performances of the two methods is the evaluation of the generalized Lomb-Scargle (GLS) periodograms (e.g., Lomb 1976;Scargle 1982;Zechmeister & Kürster 2009) on both the ∆RV * oc and ∆RV * bp residuals. The GLS periodograms for both methods are presented in Fig. 4. Thanks to the bp method, the majority of the peaks caused by stellar activity are successfully removed from the GLS periodogram (dash-dotted blue line in Fig. 4).
Conversely, the oc method is not able to properly tackle the peaks caused by active regions (dashed red line in Fig. 4), especially for periods longer than 100 days.
We used the Cramér-von-Mises (CvM) distance minimization criterion (Cramér 1928), which defines the measure we refer to as the critical value (cv), to further compare the GLS periodograms obtained from the two different methods. Going into further detail, after the GLS is obtained, the main goal is to check whether any periods are significant. A period is considered significant if its GLS periodogram peak statistically differs from the distribution that would result from the absence of periodic fluctuations (null hypothesis). To determine the significance level of the peak, this distribution needs to be known or estimated. A Beta distribution B(α, β) whose parameters are α = m−1 2 and β = n−m 2 , where n is the number of data points and m is the number of GLS parameters, is usually assumed (see e.g., Schwarzenberg-Czerny 1998;Seber & Lee 2003;Gupta & Nadarajah 2004). However, adopting those α and β values is not flexible enough, as some period values identified as significant may turn out to be false positives (see e.g., Thieler 2014;Thieler et al. 2016). In order to avoid the selection of those false positive periods, Thieler (2014) and Thieler et al. (2016) propose relaxing the assumption on the Beta distribution, by not defining its parameters "a priori", which makes the distribution more flexible. Once the assumptions on the Beta distribution are relaxed, Thieler (2014) and Thieler et al. (2016) suggest estimating the parameters of the B(α, β) by minimizing the following CvM distance: where θ is the parameter space (θ = (α, β) in our case), u = u (1) , . . . , u (n) is the vector of the ordered set of the observations (data points), F n (u) is the empirical distribution function (i.e., the empirical cumulative distribution function built from the data), and F θ (u) is the theoretical distribution function (the Beta distribution in our case). Once an estimate for θ is found (i.e., θ = (α,β)), the cv is calculated as the q √ 0.95 -quantile of the CvM-fitted Beta distribution B(α,β), where q indicates the number of period grid-points used to build the GLS periodogram. According to this criterion, the GLS peaks above the computed cv value are significant and deserve further investigation. In this work, we used the cv rather than the well known false alarm probability (FAP) because the cv is more robust (e.g., Thieler 2014; Thieler et al. 2016).
The cv computation for ∆RV * oc and ∆RV * bp yields cv oc = 0.43 and cv bp = 0.022. Peaks in the GLS below the cv are to be considered as caused by stellar activity, while more detailed analyses might be needed for peaks above the cv to discover their nature. All in all, the high cv oc value suggests quite a noisy GLS, which prevents an effective detection of exoplanetary signals. Instead, the much lower cv bp indicates a cleaner GLS, meaning that it is worth investigating even low-powered peaks, which in principle might reveal the presence of small exoplanets. We specifically checked the behavior of the GLS periodogram, by inspecting the period interval that includes the rotation period of α Cen B. Looking at Fig. 5, we note that the oc method is still unable to remove the peak at P ∼ 39 days, which is close to P rot of α Cen B (e.g., DeWarf et al. 2010). In fact, the 39-day peaks in the original and in the oc-corrected GLS periodograms both have normalized powers P RV ≈ P oc ≈ 0.17. Instead, in the bp method, the normalized power of this peak is sensibly reduced to P bp ≈ 0.027. Both P oc and P bp cannot bring us to immediately postulate the existence of a candidate exoplanet with an orbital period of 39 days. In fact P oc < cv oc = 0.43, while P bp is quite close to cv bp = 0.022, although greater. Regardless, as our main goal is cleaning the original RV time series from stellar activity, and P oc = 0.167 in the original and in the oc-corrected GLS, respectively. The normalized power of this peak is reduced at P bp = 0.027 for the bp-corrected GLS. The cv value for ∆RV * bp is displayed as the horizontal blue line and is equal to 0.022; we recall that the cv for ∆RV * oc is 0.43. it is better to deal with the bp method because it recognizes and reduces the GLS-P(P rot ) (although the corrected peak is higher than cv bp , which will encourage further investigatations to better understand its nature), rather than dealing with the oc method where the unchanged GLS-P(P rot ) (which is below cv oc ) may lead to a false negative case. As such, since P bp > cv bp , we refer the reader to Table 6 and to our discussion in Sect. 4, where we investigate the significance of all the ∆RV * bp -peaks whose normalized powers are greater than the cv threshold. We anticipate that no Keplerian-like signals have been detected.
Finally, to test the strength of the bp method, we also applied it to the RV data set analyzed by Dumusque et al. (2012) and A127, page 9 of 29 A&A 664, A127 (2022) Notes. The second to last row lists the p-values of the Wilcoxon tests, which were used to compare the rms of ∆RV * bp forD = 4 (the optimal solution) with the rms of ∆RV * bp for the other choices of D. Finally, the last row reports the BIC differences of the various D-models with respect to the optimalD-model. then by Rajpaul et al. (2015). As already noted, that data set made of 459 CCFs is a subset of the RV time series considered in this work. The GP framework presented in Rajpaul et al. (2015) contains significantly fewer free parameters than the model by Dumusque et al. (2012) (14 vs. 23). Using our bp method, we found 3 piecewise stationary segments. This makes our model comparable to the GP model of Rajpaul et al. (2015) in terms of free parameters (12 in our case), with the advantage that the linear correction we propose in Eq. (1) is simpler than any GP framework. According to our analyses, whose consequent GLS periodograms are displayed in Fig. 6, there is not any statistical evidence showing the presence of an exoplanet having an orbital period of ∼3.2 days, disproving the discovery by Dumusque et al. (2012) and confirming the conclusions presented in Rajpaul et al. (2015). In fact, our GLS normalized powers are well below the respective cv bp and cv oc thresholds for a wide neighborhood of 3 days (see Fig. 6). Therefore, we conclude that the GLS signal at ∼3.2 days visible in Fig. 4 of Dumusque et al. (2012) and announced as an exoplanet is instead likely caused by an overfitting issue.

Comparison between the optimal and the suboptimal solutions of the bp method
The bp method does not only return the bestD and the corresponding change point locations in the time series; the bp method also returns the best partition for each D = 1, . . . , D max . We can interpret those solutions as the best segmentations for each fixed D. When working with α Cen B, we foundD = 4 as the best number of breakpoints, by applying the RSS criterion to Eq. (4), combined with the BIC-based penalty function of Eq. (5). Indeed theD = 4-model has the lowest BIC among all the other models, as emphasized by the ∆BIC values reported in Table 7. We decided to further compare the results obtained from the best partition (i.e., the optimal solutionD = 4) with the other suboptimal solutions derived from the other D values. From a statistical standpoint, we found that the solutions for D = 2, D = 3 or D = 5 did not differ noticeably from the optimal one in terms of rms computed on ∆RV * bp ; the results are presented in Table 7. In particular, we used the Wilcoxon test (e.g., Zimmerman & Zumbo 1993) to compare the rms of ∆RV * bp for D = 4 with the rms of ∆RV * bp for the other available cases (null hypothesis: the compared samples are not statistically different). In detail, the periodograms have been computed from the original RV time series (solid black line), ∆RV * oc (dashed red line), and ∆RV * bp (dashdotted blue line). The plot focuses on the neighborhood of 3.2357 days (highlighted as a vertical gray line), which corresponds to the orbital period of the exoplanet announced by Dumusque et al. (2012) and then disproved by Rajpaul et al. (2015). We note that all the GLS peaks are below the cv thresholds for both of the methods.
for D = 2, 3, 4, and 5 are not statistically different for the commonly employed significance levels: α fixed = (0.01, 0.05, 0.1). However, the solution forD = 4 leads not only to the smallest rms of the ∆RV * bp time series, but also to the smallest cv among all the considered solutions, as shown in Fig. 7. We further investigated the robustness of the estimated cv by performing a bootstrap simulation study (e.g., Efron & Tibshirani 1993). The results confirmed that the cv found when usingD = 4 is significantly smaller than the cv retrieved by using suboptimal solutions; in fact the standard deviation of each cv is of the order of 10 −3 . Since all the GLS peaks below the cv threshold are ignored, dealing with low cv values means not cancelling out low GLS peaks, which might represent putative low-mass planets. As a consequence, from an astronomical standpoint, the optimal solution retrieved by the bp method is always preferable, especially when seeking super-Earths or Earth-like exoplanets.  Fig. 4 as the dash-dotted blue line). The GLS periodograms derived from the ∆RV * bp when using D = 2, 3, and 5 are displayed as a dotted blue line, a dash-dotted orange line and a dashed red line, respectively. The x-axis uses the base-10 logarithmic transformation of the period in order to improve the readability of the plot. The cv values for ∆RV * bp are displayed as black, blue, orange and red lines for the cases D = 4, 2, 3, and 5, respectively. Although the rms of ∆RV * bp for the considered cases are not statistically different, the bp-optimal solution (D = 4) leads to the smallest cv.

Simulation study based on α Cen B data
After the results were obtained by analysing real RV data of α Cen B, we decided to perform a simulation study to test how well the bp method works in detecting exoplanets. As already stressed, the RV time series of α Cen B is essentially made of stellar activity signals since no exoplanets have been detected (e.g., Rajpaul et al. 2015). Therefore, we added artificial Keplerian signals (RV K ) that would be generated by exoplanets in circular orbits to the original RV time series, we applied both the bp and oc methods to clean the data set for stellar activity (obtaining ∆RV p bp and ∆RV p oc , respectively), and we checked to which extent the two methods are able to find the simulated planets, distinguishing the planetary signals from stellar activity. Any Keplerian signal added to the RV time series translates in a shifting of the CCFs that changes according to the period, the amplitude, and the phase that is imposed. In particular, we used the following grid of values: periods from 1 to 500 days, by steps of 1 day for periods shorter than 50 days, and then by steps of 10 days for longer orbital periods, semiamplitudes from 0.1 to 15 m s −1 by steps of 0.1 m s −1 and 51 evenly sampled phases between 0 and 2π, for a total of 726 750 simulated planets.
Similar simulation studies have already been proposed by Howard et al. (2010); Mayor et al. (2011) and Simola et al. (2019). In this work, the detection of an exoplanet was tested by inspecting the GLS produced from the ∆RV p bp and ∆RV p oc data. This was implemented by first flagging as a positive finding those GLS features that are above the cv previously calculated for ∆RV * oc (cv oc = 0.43) and for ∆RV * bp (cv bp = 0.022) 4 . Then, we checked that the normalized power at P p (i.e., P p = P(P p )) was statistically comparable to the theoretical normalized powerP p . This would appear from the time series only made of the Keplerian signal of the synthetic planet that had been perturbed with a normal distribution having a mean equal to 0 and a standard deviation computed by subtracting in quadrature the rms induced by RV K from the rms of the ∆RV p time series. Because the theoretical normalized powers produced by any given exoplanetP p follow a normal distribution having constant standard deviation equal to σP = 0.028 (which can be interpreted as the inner variability of the GLS periodogram), P p is compared toP p with the commonly used Wald test (e.g., Gourieroux et al. 1982). In particular, for a given exoplanet, the test checks whether P p belongs to the 99% confidence interval inferred from N(P p , σP). This test enables us to establish whether the expected and observed RV semiamplitudes are consistent, so to declare the synthetic exoplanet as detected. This procedure is repeated for each of the 726 750 simulated planets. In order to quantify the K detection threshold, we search for the minimum RV amplitude at which, for a given orbital period, at least 90% of the planets having different phases are detected. The results of the simulation study, presented in Fig. 8, confirm the intuition that properly dividing the RV data into segments (where each segment is piecewise stationary) significantly improves the chances of detecting exoplanets that have smaller amplitudes with respect to the oc method. In fact, by applying the correction for stellar activity based on the bp method, the detection threshold is on average lower by 74% with respect to the detection threshold estimated using the oc method. The median threshold for the bp method is 2 m s −1 , while the median threshold for the oc method is 7.55 m s −1 . and therefore decided to use the already computed cv oc = 0.43 and cv bp = 0.022 throughout the following analysis. It is the minimum semiamplitude K at which a synthetic exoplanet of a given period is recovered from the two different GLS periodograms. The minimum K values are inferred from the ∆RV p bp time series (bp method, dashed blue line) and from the ∆RV p oc time series (oc method, solid red line). Blue crosses (bp method) or red crosses (oc method) highlight the null detection of any exoplanet having K ≤ 15 m s −1 . For both methods, a null detection occurs at ∼365 days, where the GLS normalized power drops. Overall, by correcting for stellar activity on each of the five temporal segments found by the bp method, the detection threshold of an exoplanet lowers by 74% on average. The two boxplots represented in the top right inset show the distributions of the minimum K for the two methods. The two distributions are statistically different as quantified by the Wilcoxon test through its p-value = 2.22 × 10 −16 . The detection threshold of an exoplanet lowers by 78% on average, when focusing on planets up to an orbital period of 250 days.
A statistical test to compare the two vectors of minimum RV amplitudes was also carried out. Assuming as null hypothesis that the two groups are not statistically different, the Wilcoxon test estimated a p-value equal to 2.22 × 10 −16 , which means that the null hypothesis is strongly rejected and that the two groups are statistically different. The detection threshold of an exoplanet lowers by 78% on average when focusing on planets up to an orbital period of 250 days. Figure 8 displays the orbital periods for which the oc method and the bp method were not able to detect any planet having K lower than our grid upper limit. In order to explain this situation, we further explored the behavior of the GLS, highlighting our findings in Sect. 3.4.2.

Period [days]
Normalized Power Fig. 9. Relation between the peak normalized power P of the GLS periodogram obtained from a synthetic Keplerian RV signal and the period of that signal (i.e., the orbital period of the synthetic exoplanet). The different colors represent the different K values imposed on the Keplerian signals, which increase from bottom to top. The ruggedness of the lines is due to the white noise that was also injected into the RV data series. Regardless of the considered K (i.e., the color), P decreases at ∼365 days (close to the revolution period of the Earth), as expected in ground-based observations.

Behavior of the GLS periodogram
We investigated the behavior of the GLS periodogram in order to understand the inability to detect exoplanets. In particular, we were interested in understanding if the unavailability of detections at a given P orb were caused by the criteria we employed for our simulation study or rather by an issue related to the GLS periodogram. We generated several Keplerian signals sampled upon the epochs of observations of α Cen B using the following grid of values 5 : periods P orb from 1 to 500 days, by steps of 1 day, semiamplitudes K from 0.1 to 15 m s −1 by steps of 0.1 m s −1 and a phase value, randomly sampled from the [0, 2π] interval. Each vector of RV was perturbed by using a normal distribution having a mean equal to 0 and a standard deviation of 1 m s −1 , which is the error expected by HARPS. At this stage, as vector of RV, just the Keplerian signal and the error term are used.
When plotting P as a function of P orb for any given value of K, rather than finding constant behavior, there is a problematic behavior when considering short orbital periods. Moreover, we see a systematic decrease of P at some specific P orb values, which is caused by the data sampling due to α Cen B visibility. In particular, a clear decrease occurs at P orb ∼365 days, as expected in ground-based observations. These conclusions are confirmed in our simulation study that aims to establish the detection threshold of the bp and oc methods, as both the methods have some issues in detecting planets that have a short orbital period, and both methods are unable to detect planets with P orb ∼365 days (Fig. 9).

Comparison of the oc and the bp methods applied to
HD 215152, HD 10700, and HD 192310 As done for the RV time series of α Cen B, we repeated the comparison between the oc and the bp methods for three other main sequence stars: HD 215152 (K3V, Delisle et al. 2018), HD 10700 (G8V, Feng et al. 2017b), andHD 192310 (K2V, Pepe et al. 2011). Unlike α Cen B, which does not host any known planet, HD 215152 and HD 10700 host four planets (Delisle et al. 2018;Feng et al. 2017b, each), while HD 192310 hosts two planets (Pepe et al. 2011). Therefore, to obtain ∆RV * -like indicators that are consistent with those derived for α Cen B, we should also remove the exoplanetary signals from these RV time series, after applying the correction for stellar activity. Since the discovered exoplanets orbiting HD 215152 and HD 10700 are at the level of instrumental precision, removing those planetary signals would have led to biased results. For this reason, in the following analyses we decided to remove the planetary signals only for the HD 192310 case. We used the adaptive Markov chain Monte-Carlo (A-MCMC) algorithm (Haario et al. 2005) to estimate the orbital parameters of the two known exoplanets orbiting HD 192310. In detail, following a step-by-step procedure, first we corrected the original RV time series by using both the oc and the bp methods. As a result, a clear peak due to the exoplanet with P orb ≈ 75 days appeared in the two GLS periodograms (although only in case of the bp method P bp (75 d) > cv bp ; top panel of Fig. 10). Therefore, we estimated the orbital parameters of this exoplanet by using the A-MCMC algorithm, and removed its signal (RV P:75 ) from the ∆RV * time series. After producing new GLS periodograms from the ∆RV * − RV P:75 time series, for both methods we saw another peak above the respective cv values, which was compatible with P orb ≈ 526 days of the second exoplanet (second panel of Fig. 10). Then we used the A-MCMC algorithm again to estimate the orbital parameters of the exoplanet with P orb ≈ 526 days, we removed its signal (RV P:526 ) from the previously used RV time series, and we produced new GLS periodograms based on the ∆RV * − RV P:75 − RV P:526 time series. The inspection of these GLS periodograms did not show any other significant peak (third panel of Fig. 10) and we concluded that only the bp method finds both the two known exoplanets following the cvthreshold criterion. Finally, we launched a global A-MCMC run  Notes. The relative uncertainties are specified within brackets, while the ∆ σ column refers to the statistical tension between the generic pair (x 1 ± σ 1 ; accounting for both planets at the same time to obtain unbiased exoplanetary parameters, which were used to build our final GLS periodograms, where both the planetary signals were subtracted from ∆RV * (bottom panel of Fig. 10).
Our outcomes for the HD 192310 planetary system are listed in Table 8 and compared with the results obtained by Pepe et al. (2011). The precision we got for the planet b parameters is comparable with that obtained by Pepe et al. (2011) and all those parameters are consistent within ∼1.5σ, except for the orbital period P b . Our P b = 74.06 +0.10 −0.09 days is 4.8σ away from the Pepe et al. (2011) estimate. However, it is in excellent agreement with the recent revision by Rosenthal et al. (2021), who found a period of 74.062 ± 0.085 days. Regarding planet c, instead, we sensibly improve the precision on the determined parameters (from a factor of two up to a factor of four) with respect to the results provided by Pepe et al. (2011). All the respective estimates are within 2σ, except for K c , for which we infer a higher value after applying our bp method for cleaning the RV time series. Even if a detailed analysis of the orbital parameters is beyond the scope of this paper, we note that our estimates for HD 192310 c are the most precise available in the literature so far 6 . Coming to the effectiveness of modeling the stellar activity, when compared to the rms of the ∆RV * oc time series, the rms of ∆RV * bp is lower by 0.21 m s −1 (HD 215152), 0.06 m s −1 (HD 10700), and 0.1 m s −1 (HD 192310, after the planets removal). This means that using Eq. (6) rather than Eq. (7) for correcting RV, we better explain the RV-variability by gaining 9.5%, 4.3%, and 6.6% for HD 215152, HD 10700, and HD 192310, respectively. The bp method has to be preferred also following the BIC minimization criterion for model selection; in fact, ∆BIC = BIC oc − BIC bp = +79, +503, and +174 for HD 215152, HD 10700, and HD 192310, respectively. The number of available data points per time series is significantly smaller if compared to the α Cen B time series, especially for HD 215152 and HD 192310. It seems therefore reasonable that the bp method finds a lower number of piecewise-stationary segments as optimal solutions for those stars (D = 1 and 2 for HD 215152 and HD 192310, respectively). Instead, concerning HD 10700, a larger number of observations spread on a longer baseline is available, which yields toD = 5 optimal breakpoints. 6 According to the Nasa Exoplanet Archive (https:// exoplanetarchive.ipac.caltech.edu/overview/HD192310), where only the set of parameters by Pepe et al. (2011) is provided.
As a final note, since HD 10700 is a quiet star, the bp vs. oc improvement in modeling the stellar activity is lower (i.e., we registered a lower gain when comparing the ∆RV * bp vs. ∆RV * bp rms, as expected). However, the bp method still models the stellar activity better than the oc method and the ∆BIC confirms a strong preference for the bp method. Overall, the obtained results confirm our assumption that the CPD algorithm is particularly effective when focusing on active stars. Relevant tables and plots synthesizing our results can be found in Appendix A.

Discussion
We tested the bp method by using real measurements taken from four stars, carrying out comparisons with the oc method that considers a single correction for stellar activity on the entire time series. Before performing the comparisons, the data were properly cleaned by removing the outliers from the set of RV values. The results suggest that properly dividing the RV time series into segments (where each segment is piecewise stationary) is a helpful operation when trying to account for higher variations in RV caused by stellar activity. For all the considered stars, the rms on ∆RV * bp are smaller than the rms on the corresponding ∆RV * oc . The stronger the activity signals within the time series, the more effective the bp method is at modeling the RV variations. In addition, the longer the time series, the more likely sensible variations in stellar activity have occurred, and so the bp method is particularly suitable to detrend RV data. This is especially the case of α Cen B, which is characterized by strong solar-like activity signals (e.g., Thompson et al. 2017;Dumusque 2018) . By dividing the RV measurements of α Cen B into five piecewise stationary segments, we reach an rms of ∆RV * bp of 1.75 m s −1 (to be compared with an rms of ∆RV * oc of 3.48 m s −1 resulting from the application of the oc method). Also, the activity indicators' variability within each segment is smaller than the variability between each segment, as displayed in Fig. 1 and summarized by Table 3, which provides strong evidence that the segmentation proposed by the bp method captures those time series locations where the stellar activity changes significantly. As a consequence, the bp method is able to interpret a larger fraction of RV variability in terms of stellar activity: 89% vs. 40% for the oc method.
The α Cen B GLS periodogram obtained with the bp method shows a much smaller number of peaks caused by active regions A127, page 15 of 29 A&A 664, A127 (2022) with respect to the GLS periodogram obtained by performing the oc method. When producing the GLS periodogram from the ∆RV * oc time series, we found that there were no peaks above cv oc . Conversely, when considering the ∆RV * bp time series, 11 peaks in the GLS were above the critical value of cv bp (see Table 6). We further investigated the nature of those peaks to check whether they could be produced by exoplanets. For a given exoplanet candidate, starting from ∆RV * bp and assuming a Keplerian model, we used the A-MCMC algorithm to retrieve the posterior distribution of the parameters of interest. Then, we used the marginal posterior means for each parameter in order to estimate the RVvariations caused by the candidate planet. Finally, we compared the rms of ∆RV * bp with the new rms obtained by subtracting to ∆RV * bp the RV Keplerian signal that would be caused by the candidate planet. We did not find any strong signal suggesting the possible presence of an exoplanet, confirming the conclusions provided in Rajpaul et al. (2015).
When repeating the comparative analyses for HD 215152, HD 10700, and HD 192310, the bp method produces ∆RV * bp time series having rms lower by 0.21 m s −1 for HD 215152, 0.06 m s −1 for HD 10700, and 0.1 m s −1 for HD 192310 when compared to the rms of ∆RV * oc . In other words, the bp method increases the fraction of RV variability that can be explained in terms of stellar activity by 9.5%, 4.3%, and 6.6% for HD 215152, HD 10700, and HD 192310, respectively. The improvement given by the bp method is less evident in these three cases, as these stars are less active than α Cen B and their time series span a shorter temporal range. Since these time series already appear to be piecewise stationary as a whole, the bp and oc methods are comparable. However, thanks to its better cleaning performance, only the bp method was able to detect both the exoplanets hosted by HD 192310, following our cv-threshold criterion. In particular, we essentially confirm the orbital parameters estimates already available in the literature, but we sensibly improve the precision in the case of HD 192310 c. Given the key role of the planetary mass when studying the composition of exoplanets, we emphasize that our M p sin i estimates are affected by an error of ∼5% for both planet b and planet c. We reached the same precision level of Pepe et al. (2011) for planet b, while we improved the precision of M p sin i of HD 192310 c by a factor of approximately four.
To further evaluate the ability of the bp and oc methods to detect exoplanets, we designed a simulation study starting from α Cen B data. We added several synthetic Keplerian signals to the time series and checked the exoplanet detection effectiveness of both the bp and oc methods when dealing with RV data points contaminated by solar-like activity signals. We found that the bp method lowers the detection threshold (i.e., the minimum K of the Keplerian signal at which a planet is detected from the inspection of the GLS periodogram) with respect to the oc method by 74% when considering planets up to an orbital period of 500 days.

Conclusion
Stellar activity, in the form of active regions evolving on a star's photosphere, has so far been the major obstacle for the detection and the characterization of Earth-like exoplanets when using the RV method. Spots and faculae cause variations in the shape and in the width of the CCF, changing the correlations between RV and the indicators of stellar activity, such as A, γ, and FWHM SN . Since an exoplanet would not change the shape of the CCF, but just its barycenter, a common strategy to account for stellar activity is to employ a linear correction of the RV time series involving the activity parameters. In fact, it is well known that variations in the correlations between RV and these activity indicators suggest the presence of active regions evolving on the stellar photosphere over time. A simple way to model the changes in RV caused by stellar activity is provided by Eq. (1). Since RV surveys often spread over years of measurements, it seems reasonable to assume that the stellar activity level changes multiple times during the observational period. Rather than using an overall correction for stellar activity on the entire time series, we still rely on Eq. (1), but we propose performing multiple corrections by suitably dividing the overall time series into segments. The number of segments depends on how often the correlations between RV and the activity indicators significantly change.
In order to estimate the time series locations where the dependence of stellar activity upon activity indicators significantly changes, we draw attention to the family of CPD algorithms.
In particular, in this paper we used the CPD-based bp method (e.g., Bai & Perron 2003) to properly model the variations of RV caused by stellar activity. We compared the effectiveness of the bp method with the commonly employed oc method, by using real observations taken on four different stars. The results show that identifying the locations in the RV data where the correlations between RV and the indicators of stellar activity significantly change, produces much cleaner RV time series as we model the stellar activity signals on each of the piecewise stationary segments. The GLS periodograms are then less contaminated by the presence of spurious periodical peaks caused by stellar activity. As demonstrated by our simulation study on α Cen B, the bp method was able to detect exoplanets that produce RV amplitudes 74% smaller than those detected by the oc method.
Finally, we note that the bp method is most effective when working with active stars whose RV time series are made of several hundreds of data points. In fact, the longer the time series, the more likely sensible variations in stellar activity have occurred, suggesting that the bp technique is a suitable statistical tool for removing activity-induced variations from the RV data.         Fig. 3, but for HD 192310. We recall that the two planetary signals were removed from our analyses after the bp method and the oc method were used, in order to better focus on the residuals of the two methods.
A127, page 28 of 29 Table A.10. Normalized power P and cv values inferred from the bp-corrected GLS (P bp , cv bp ) and oc-corrected GLS (P oc , cv oc ) for HD 192310. They are evaluated at the stellar rotation period (P rot ) and at the at the orbital periods (P orb ) of the hosted exoplanets found by Pepe et al. (2011). The first two rows show P inferred from the ∆RV * signals, while in the last row P are computed after removing the Keplerian signal of P orb,1 = 75 days from the ∆RV * signals. Both the two planets would be visible in the GLS after applying the bp method and the iterative removal of planetary signals, as P bp are always above the respective cv bp values.