© ESO 2019

1 Introduction

When working with radial-velocity data (RVs), one of the main limitations to the detection of low-mass exoplanets is no longer the precision of the instruments used, but the different sources of variability induced by the stars (e.g., Feng et al. 2017a; Dumusque et al. 2017; Rajpaul et al. 2015; Robertson et al. 2014). Stellar oscillations, granulation phenomena, and stellar activity can all induce apparent RV signals that are above the meter-per-second (m s⁻¹) precision (e.g., Saar & Donahue 1997; Queloz et al. 2001; Desort et al. 2007; Dumusque et al. 2011; Dumusque 2016) reached by the best high-resolution spectrographs (HARPS, HARPS-N; Mayor et al. 2003; Cosentino et al. 2012). It is therefore mandatory to better understand stellar signals and to develop methods to correct for these signals, if in the near future we want to detect or confirm an Earth-twin planet using the RV technique. This is even more true now that instruments like the Echelle SPectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2014) and the EXtreme PREcision Spectrometer (EXPRES; Fischer et al. 2016) should reach the precision and stability to detect such signals. However, if solutions are not found to mitigate the impact of stellar activity, the detection or confirmation of potential Earth-twins will be extremely challenging and false detections could plague the field.

One of the most challenging stellar signals to characterize and correct for are the signals induced by stellar activity. Stellar activity is responsible for creating magnetic regions on the surface of stars, and those regions change locally the temperature and convection, which can induce spurious RVs variations (e.g., Meunier et al. 2010; Dumusque et al. 2014; Borgniet et al. 2015). In theory, it should be easy to differentiate between the Doppler-shift induced by a planet, which shifts the entire stellar spectrum, and the effect of stellar activity, which modifies the shape of spectral lines and by doing so creates a spurious shift of the stellar spectrum (Saar & Donahue 1997; Hatzes 2002; Kurster et al. 2003; Lindegren & Dravins 2003; Desort et al. 2007; Lagrange et al. 2010; Meunier et al. 2010; Dumusque et al. 2014). However, on quiet GKM dwarfs, the main targets for precise RVs measurements, stellar activity can induce signals of a few m s⁻¹. This corresponds physically to variations smaller than 1/100th of a pixel on the detector, making changes in the shape of the spectral lines challenging to detect. In order to measure such tiny variations, a common approach is to average the information of all the lines in the spectrum by cross correlating the stellar spectrum with a synthetic or an observed stellar template (Baranne et al. 1996; Pepe et al. 2002; Anglada-Escudé & Butler 2012). The result of this operation gives us the cross-correlation function (CCF). To measure the Doppler-shift between different spectra, and therefore to retrieve the RVs of a star as a function of time, the variations of the CCF barycenter are calculated. The barycenter is generally estimated by the mean of a Normal density shape fit to the CCF. Variations in line shape between different spectra, which indicate the presence of signals induced by stellar activity, are measured by analyzing different parameters of models fit to the CCF. Usually the width of theCCF is estimated using the full width half maximum (FWHM) of the fitted Normal density and its asymmetry by calculating the CCF bisector and measuring the bisector inverse slope span (BIS SPAN; Queloz et al. 2001).

If a spurious RV signal is induced by activity, generally a strong correlation is observed between the RV and chromospheric activity indicators like log () or H-α (Boisse et al. 2009; Dumusque et al. 2012; Robertson et al. 2014), but also between the RV and FWHM of the CCF or its BIS SPAN (Queloz et al. 2001, 2009; Boisse et al. 2009; Dumusque 2016). Therefore, a common strategy when fitting a Keplerian signal to a set of estimated RVs when searching for a planet is to include linear terms in the model to account for activity, such as the log (), FWHM, and BIS SPAN (Dumusque et al. 2017; Feng et al. 2017a). It is also common to add a Gaussian process to the model to account for the correlated noise induced by stellar activity. The hyperparameters of the Gaussian process can be trained on different activity indicators (Haywood et al. 2014; Rajpaul et al. 2015) or directly on the RVs (Faria et al. 2016). It is therefore essential for mitigating stellar activity to obtain activity indicators that are the most correlated with the RVs but also for which we can obtain the best precision.

Several indicators have been developed that can be more sensitive to line asymmetry than the BIS SPAN. In Boisse et al. (2011), the authors developed V_span, which is thedifference between the RV measured by fitting a Normal density to the upper and lower parts of the CCF.This CCF asymmetry parameter is shown to be more sensitive than the BIS SPAN at a low signal-to-noise ratio (S/N). Figueira et al. (2013) studied the use of new indicators, BIS-, BIS+, bi-Gauss, and V_asy. The authorswere able to show that when using bi-Gauss, the amplitude in asymmetry is 30% larger than when using BIS SPAN, therefore allowing the detection of lower levels of activity. They also demonstrated that V_asy seems to be a better indicator of line asymmetry at high S/N, as its correlation with RV is significantly stronger than any other correlation between the previously proposed asymmetry indicators and RV.

In all the methods described above, except the bi-Gauss, RV, and FWHM are derived using a Normal density fitted to the CCF, and the asymmetry is estimated using a separate approach. In this paper we propose to use a Skew Normal (SN) density to estimate with a single fit of the CCF the RV, FWHM, and asymmetry of the CCF, as this function includes a Skewness parameter (Azzalini 1985).

The paper is organized as follows. In Sect. 2 we introduce the SN density, describe its applicability for modeling the CCF, and study how the SN parameters relate to the mean of the Normal density, FWHM, and BIS SPAN of the CCF. In Sect. 3 we propose a linear model to correct for stellar activity signals in RVs, which extends the linear models previously proposed for this purpose (e.g., Dumusque et al. 2017; Feng et al. 2017a). In Sect. 4 the performance of the SN fit to the CCF is investigated using simulations coming from the Spot Oscillation And Planet 2.0 code (SOAP 2.0; Dumusque et al. 2014), followed by an analysis of real observations in Sect. 5. In Sect. 6 error bars are computed for the different estimated CCF parameters. Finally, the discussion of the results and conclusions are presented in Sects. 7 and 8, respectively.

2 Skew Normal distribution

The SN distribution is a class of probability distributions that includes the Normal distribution as a special case (Azzalini 1985). The SN distribution has, in addition to a location and a scale parameter analogous to the mean and standard deviation of the Normal distribution, a third parameter that describes the Skewness (i.e., the asymmetry) of the distribution. Considering a random variable (where is the real line), which follows a SN distribution with location parameter , scale parameter (i.e., the positive real line), and Skewness parameter , its density at some value y ∈ Y can be written as (1)

where ϕ and Φ are, respectively, the density function and distribution function of a standard Normal distribution¹. The Skewness parameter α quantifies the asymmetry of the SN. Examples of SN densities under different Skewness parameter values and the same location and scale parameters (ξ = 0 and ω = 1) are shown in Fig. 1. A usual Normal distribution is the special case of the SN distribution when the Skewness parameter α is equal to zero².

For reasons related to the interpretation of the parameters in Eq. (1) and computational issues with estimating α near 0, a different parametrization is used in this work, which is referred to as the centered parametrization (CP). This CP is much closer to the parametrization of a Normal distribution, as it uses a mean parameter μ, a variance parameter σ², and a Skewness parameter γ. In order to define the CP, we need to express the CP parameters (μ, σ², γ) as a function of (ξ, ω², α). This can be done using the following relations: (2)

where (e.g., Arellano-Valle & Azzalini 2008).

By using Eq. (2), the new set of parameters (μ, σ², γ) provides a clearer interpretation of the behavior of the SN distribution. For the α values used in Fig. 1, the corresponding values of (μ, σ², γ) are shown in Table 1. In particular, μ and σ² are the actual mean and variance of the distribution, rather than simply a location and scale parameter, and γ provides a measure of the Skewness of the SN. Along with the mean of the SN, we consider the median of the distribution as a measure of its barycenter. See Table 1 for the medians of the SN densities shown in Fig. 1.

Further details about the parametrization from Eq. (1), called the direct parametrization (DP), the CP, and general statistical properties of the SN can be found in Azzalini & Capitanio (2014).

To fit the CCF using a SN density shape, we use a least-squares algorithm and the following model: (3)

where Cis an unknown offset for the continuum of the CCF, A is the unknown amplitude of the CCF, sometimes referred to as the CCF contrast, and μ, σ², and γ are the mean, variance, and Skewness of the SN, respectively, as defined above. The values x₁, …, x_n are the different values of the x-axis of the CCF, generally in velocity units (e.g., m s⁻¹).

When using a Normal density shape model for the CCF, the estimated mean is used as the estimated RV and the FWHM³ is used to quantify the width of the CCF. Because the Normal density is symmetric, the Skewness is not defined and therefore a separate approach is necessary to estimate the Skewness of the CCF. An estimated Skewness parameter is generally obtained by calculating the BIS SPAN of the CCF (see Sect. 1 and, e.g., Queloz et al. 2001).

With the proposed SN approach, we propose two estimators of the RV: the mean and median of the SN model fit (referred to as SN mean RV and SN median RV, respectively), and present advantages and limitations for both of these choices in Sects. 5 and 6. The width of the SN, SN FWHM is defined in the same way as for the Normal density⁴, and finally the Skewness of the CCF is estimated by the γ parameter.

To evaluate the strength of the correlation between the estimated RVs and the different stellar activity indicators, we calculated the Pearson correlation coefficient, R, which in its general form is defined as (4)

where x and y are two quantitative variables, cov(x, y) indicates the covariance between x and y, and σ(x) and σ(y) represent their standard deviations. A p-value for the statistical test with null hypothesis H₀: R = 0 is also generally provided.

	Fig. 1 Density function of a random variable Y following a SN distribution SN(ξ, ω, α) with location parameter ξ = 0, scale parameter ω = 1 and differentvalues of the Skewness parameter α indicated by different colors and line types. We note that the solid black line has an α = 0 making it a Normal distribution.
Open with DEXTER

3 Radial-velocity correction for stellar activity

Exoplanets produce a Doppler shift of the entire stellar spectrum. On the contrary stellar activity and, in particular, the presence of active regions on the stellar photosphere, do not produce blueshifts or redshifts of the entire stellar spectrum but can create spurious RVsignals by modifying the shape of spectral lines. To track these variations in the shape of the spectral lines, a common approach consists in using the FWHM, BIS SPAN, or other indicators such as those introduced in Boisse et al. (2011) or Figueira et al. (2013), which provide information on the width and asymmetry of the CCF. A strong correlation between the estimated RVs and one or more of these parameters provides an indication that stellar activity signals may be affecting the measurements.

When fitting for planetary signals in RV data, it is common to include linear dependencies with the BIS SPAN and FWHM to take into account the signal induced by stellar activity (e.g., Dumusque et al. 2017; Feng et al. 2017a). We propose to add additional parameters in the model to correct for stellar activity: first, an amplitude parameter A of the CCF (referred to as the CCF contrast) and second, an interaction term for γ and SN FWHM (or the BIS SPAN and the FWHM in the Normal case). The stellar activity correction we propose can therefore be written as (5)

where β₀ is the intercept and ϵ is the random error with mean equal to 0 and covariance matrix equal to σ²I (I defined as the identity matrix). The contrast parameter A accounts for the presence of a spot on the stellar surface, which produces a change in the amplitude of the CCF, in addition to changes in asymmetry or width (see, e.g., Fig. 2 in Dumusque et al. 2014). The benefits of including a variable that quantifies the interaction between γ and SN FWHM (or BIS SPAN and FWHM) is better understood through the results of the examples presented in Sect. 4. This interaction term can account for possible interactions between SN FWHM (or FWHM) and γ (or BIS SPAN), meaning that the association of each variable with the response, RV_activity, depends also on the other variable.

The proposed model is analyzed using statistical tests on the parameters β₀, β₁, β₂, β₃, and β₄ where the null hypothesis is H₀: β_i = 0, for i = 0, …, 4. The significance level for the tests is set at 0.05. The coefficient ofdetermination, R², is used to assess how well the proposed linear combination of variables accounts for the variability of RV_activity.

The proposed function defined in Eq. (5) is the result of statistical and astronomical considerations. In particular, we evaluated the correlations between the model covariates, since high correlations can lead to a non-invertible design matrix resulting in invalid parameter estimates. This problem is known in statistics as multicollinearity and more discussion can be found in Belsley (1991). In the analysis of real data presented in this work, we never observed a correlation coefficient exceeding 0.66 between theasymmetry and width parameters and therefore the problem of multicollinearity appeared tobe mitigated. We also investigated the statistical significance of the interactions between A and the width of the CCF, and A and the asymmetry of the CCF; those two interactions were not statistically significant.

4 Simulation study

In order to evaluate the performance of the proposed SN approach to modeling the CCF and the benefit of using the proposed correction for stellar activity (see Eq. (5)), we begin by considering a simulation study using spectra generated from the SOAP code (Dumusque et al. 2014).

For a given configuration of spots and faculae on the stellar surface, SOAP 2.0 outputs simulated CCFs as a function of rotational phase. The code also returns the RV and FWHM estimated using a Normal model for the CCF and the BIS SPAN obtained by calculating the bisector of the CCF.

For the simulations discussed below, a star similar to the Sun was modeled with a solar disk of one solar radius seen equator-on, and a stellar rotational period set to 25.0 days. The stellar effective temperature is set to 5778 K, and a quadratic limb-darkening relation with linear and quadratic coefficients 0.29 and 0.34, respectively (Oshagh et al. 2013; Claret & Bloemen 2011). The temperature difference of the spot with the photosphere is ΔT_spot = −663 K and the temperature difference of the facula depends on the center-to-limb angle θ, ΔT_facula = 250.9 − 407.7 cosθ + 190.9 cos²θ K (Meunier et al. 2010). In order to make the result of the simulations more comparable to real data obtained with the HARPS spectrograph discussed in Sect. 5, the SOAP 2.0 CCFs were generated with a width of 40 km s⁻¹ and considering initial spectra with resolution of Res = 115 000. For the simulations presented in this section, we considered a S/N of 100.

4.1 Facula

To see the impact of a facula on the parameters of different models of the CCF, we simulated the effect of an equatorial facula of size 3% relative to the visible stellar hemisphere. The facula is face-on when the phase is 0. We note that a 3% faculae is relatively large for the Sun; at maximum activity, large faculae generally have a size of 1% (e.g., Borgniet et al. 2015). In Fig. 2, we compare the barycentric variation of the CCF as measured when fitting a Normal density and using its mean (N mean RV), and when fitting a SN density and taking its mean (SN mean RV) or median (SN median RV). We see that all the different estimates of the CCF barycenter present a signal of similar amplitude, however the signal obtained with SN mean RV is notably different from the two others with a maximum amplitude happening at a different phase.

Correlations between the different RV estimates and the different CCF asymmetry or width estimates are shown in Fig. 3. The correlation between γ and SN mean RV, and γ and SN median RV are weaker than the correlation between BIS SPAN and RV and have Pearson correlation coefficient values of R = −0.11, −0.55, and −0.61, respectively. Therefore, itseems that when looking at the CCF asymmetry, fitting a SN density to the CCF does not help. However, when looking at the correlations between the width of the CCF and the RV, the stronger correlation can be found between SN FWHM and SN mean RV (R =0.92). Then the correlation follows between SN FWHM and SN median RV (R = 0.65) and finally FWHM and N mean RV (R = 0.47). Over all the different correlations in the case of a facule, the strongest is found between SN FWHM and SN mean RV, which demonstrates that fitting a SN density to the CCF can be helpful to better probe stellar activity.

The RVvariations shown in Fig. 2 are caused only by stellar activity, in this case a facula. We applied the activity correction proposed in Eq. (5) to evaluate the performance of the model in this setting. The results of this correction are shown in Fig. 4 and the statistical tests on the coefficients involved in Eq. (5) are summarized in Table 2. The proposed correction for stellar activity is able to account for the majority of the activity signal created by a facula, with a R² larger than 0.98. In addition, the proposed linear model allows us to reduce the activity effect from a RV rms of 8.02 to 1.02 m s⁻¹ when considering the CCF contrast, FWHM, and BIS SPAN. When using the parameters derived from the SN, the improvements for SN median RV and SN mean RV are 7.07 m s⁻¹ down to 0.88 m s⁻¹ and 5.9 m s⁻¹ down to 0.88 m s⁻¹, respectively.When comparing the correction proposed in Eq. (5) with what is generally used (i.e., a linear combination of only the asymmetry and width parameters), we see that the proposed correction is able to reduce the rms of the RV residuals by an additional 14.5–15%. Looking at the significance of the coefficients in Table 2, we observe that all the parameters corresponding to the SN variables are statistically significant. When using the Normal variables, the parameters β₃ and β₄ are not statistically helpful for the correction.

Fig. 2 RV estimates for N mean RV (red dashed line), SN mean RV (black line), and SN median RV (cyan dotted-dashed line). In this case, the CCFs were generated using SOAP 2.0 with an equatorial3% facula on the simulated Sun. The star does one full rotation between phase −0.5 and 0.5; the facula are seen face-on for phase 0. The variations observed in SN mean RV are notably different from the variations measured in SN median RV and N mean RV.

Open with DEXTER

4.2 Spot

Next we consider the effects on the CCF model parameters due to the presence of an equatorial spot of size 1% relative to the visible stellar hemisphere. The spot is face-on when the phase is 0. We note that this is a large spot for the Sun, as large spots are generally onthe order of 0.1% (e.g., Borgniet et al. 2015). In Fig. 5, we show the barycentric variations of the CCF induced by this simulated spot. In contrast to the case of the facula, all the different estimates of the CCF barycenter for the spot have the same shape in variation; however, the amplitude for SN mean RV is slightly smaller.

Figure 6 shows the correlations between the asymmetry parameters and the different estimates of the CCF barycenter (i.e., SN mean RV, SN median RV, and N mean RV). The correlation between γ and SN median RV and the correlation between the BIS SPAN and N mean RV are the strongest (R = −0.81) followed by the γ – SN mean RV correlation (R = −0.76). The correlations between the width and CCF barycenter draw circles and no significant correlation is observed. Unlike the facula scenario, when considering a spot simulated from SOAP 2.0 and a S/N of 100 the SN parameters do not appear to better probe stellar activity than the Normal parameters.

As before, the original RV estimates are corrected using Eq. (5). The results of this correction are indicated in Fig. 7 and the statistical tests on the coefficients involved in Eq. (5) are summarized in Table 3. In the case of a spot, the proposed correction is not able to perform as well as for the facula, and R² values for the linear combination are between 0.7 and 0.8. The correction is able to mitigate stellar activity from a N mean RV rms of 6.14 m s⁻¹ down to 3.04 m s⁻¹. The corresponding values for SN median RV and SN mean RV are 5.85 m s⁻¹ down to 2.74 m s⁻¹ and 5.27 m s⁻¹ down to 2.74 m s⁻¹, respectively. When comparing the activity correction proposed in this paper with what is commonly used, which means only a linear combination of the width and asymmetry of the CCF, we see that our solution is capable of reducing the RV residual rms by 5.3–5.8%. The Normal or SN parameters involved in Eq. (5) are statistically significant to explain the activity signal as seen in Table 3, except the width of the CCF β₃ and the interaction term β₄.

Fig. 3 Left panels: correlations between the different asymmetry parameters and their corresponding RV estimates in the case of an equatorial 3% facula on the simulated Sun. Right panels: correlations between the different width parameters and their corresponding RV estimates for the same facula. In the presence of a facula, both the shape and the width of the CCF change as the star rotates, producing statistically significant correlations for all the cases except for the correlation between SN mean RV and SN GAMMA (P-value = 0.27).

Open with DEXTER

4.3 Spot and planet

The final simulation includes a planetary signal influencing the CCF along with the 1% spot modeled previously (see Sect. 4.2). The purpose of this example is to check if we are able to disentangle these two different sources of variations when using the parameters derived using a Normal or a SN model for the CCF. In this scenario the planet is injected with a semi-amplitude of 10 m s⁻¹ with no eccentricity and with a period corresponding to one-third of the stellar rotational period, i.e., one-third of 25 days.

Figure 8 shows the variations observed in the CCF barycenter parameters. As in the case of the spot, all RV estimates show similar variations, and SN mean RV shows a slightly smaller amplitude.

The correlations between the different CCF parameters are represented in Fig. 9. The correlations are weaker than in the case of the spot due to the planet inducing changes in RV without affecting the width or the asymmetry of the CCF. However, the order of the strength of the correlations between the CCF asymmetry parameters and RV are comparable with those obtained for the spot-only model: γ–SN median RV has the strongest correlation (R = −0.7), followed by the correlation between BIS SPAN–N mean RV (R = −0.69), and finally by the correlation between γ–SN mean RV (R = −0.67). The patterns observed in the width-RV phase-space correlations in Fig. 9 follow a circle, similar to the spot-only model; no substantial correlation is observed between those two parameters.

In orderto correct the estimated RVs from the spurious variation caused by the spot, the proposed model for correcting the activity is added to a planetary signal model that takes into account the RV variations caused by a planet. The observed RVs can therefore be modeled as a combination of the activity and planetary signals, i.e., (6)

where RV_activity can be found in Eq. (5), and RV_planet, in the case with no eccentricity, can be modeled by the following sinusoidal function: (7)

with amplitude K, orbital period P, and an epoch at the periapsis t₀. The previous three unknown parameters define the planetary orbit.

The proposed model from Eq. (6) was fitted to the RV data and the results of the estimated model are summarized in Table 4. Except for the intercept β₀, width of the CCF β₃, and interaction term β₄ that evaluates the interaction between the width and asymmetry of the CCF, all the other Normal or SN parameters are significantly useful to explain the RV variation induced by a spot plus a planet. We also observe that the RV residuals, once corrected for stellar activity and the presence of the planet, are smaller in terms of rms when using the SN variables (rms = 2.29 m s⁻¹) rather than the Normal variables (rms = 2.80 m s⁻¹).

Fig. 4 Top panels: spurious estimated RVs (black dots) caused by a facula in the simulated data using a Normal and a SN model, the estimated RVs using Eq. (5) (red pluses), and the estimated RVs using the usual correction for stellar activity (green triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (red pluses) and the residuals from the usual correction (green triangles). The standard deviations are also reported in the legend, and the residuals have a smaller systematic component when using the proposed model of Eq. (5) compared to the usual model. The tests of statistical significance on the parameters are presented in Table 2.

Open with DEXTER

	Fig. 5 RV estimates for N mean RV (red dashed line), SN mean RV (black line) or SN median RV (cyan dotted-dashed line) using CCFs generated from SOAP 2.0 with an equatorial 1% spot on the simulated Sun. The star does one full rotation between phase −0.5 and 0.5; the spot is seen face-on at phase 0. The SN mean RV seems to have the smallest spurious variations caused by the spot.
Open with DEXTER

Fig. 6 Left panels: correlations between the different asymmetry parameters and their corresponding RV estimates in the case of an equatorial 1% spot on the simulated Sun. Right panels: correlations between the different width parameters and their corresponding RV estimates for the same spot. In the presence of a spot, both the shape and width of the CCF change as the star rotates. However, only the asymmetry produces a statistically significant correlation with the different RV estimates. The width parameters and their corresponding RV estimates present weak correlations and, in general, much weaker correlations compared to the results obtained when an equatorial 3% facula is present on the simulated Sun.

Open with DEXTER

Fig. 7 Top panels: spurious estimated RVs (black dots) caused by a spot in the simulated data, the estimated RVs using Eq. (5) (red pluses), and the estimated RVs using the usual correction for stellar activity (green triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (red pluses) and the residuals from the usual correction (green triangles). The standard deviations are also reported in the legend, and the residuals have a smaller systematic component when using the proposed model compared to the usual model. The tests of statisticalsignificance on the parameters are presented in Table 3.

Open with DEXTER

5 Real data application

The analyses of the simulated SOAP 2.0 data in the previous section were helpful in assessing the performance of the proposed methodology in a setting where the ground truth was known. We found with the simulated data that the parameters derived by the SN had correlations comparable to those derived by the Normal; in the case of a facula, the SN parameters had higher correlations for the FWHM than those of the Normal. In this section we present an analysis conducted on real observations, in particular, the star Alpha Centauri B, and compare the performance when fitting a CCF using the SN density defined in Sect. 2 with the commonly employed approach based on fitting a Normal density for estimating the RV and retrieving the BIS SPAN for evaluating the asymmetry of the CCF. Four other stars have been analyzed with the proposed method and details can be found in Appendix A. For all the stars considered in the presented work, only CCFs that were derived from spectra that had at least a S/N of 10 at 550 nm were selected.

5.1 Comparison for Alpha Centauri B of the different CCF parameters derived with Normal and Skew Normal

A total of 1808 CCFs that were derived from the spectra of Alpha Centauri B taken in 2010 by the HARPS spectrograph have been analyzed. We note that more observations were carried out that year, however only the data that were not significantly affected by contamination from Alpha Centauri A were used (see Dumusque et al. 2012). The selected observations represent arguably, among all RV data existing, the best sampled and most precise RV data set showing strong solar-like activity signals (Dumusque 2018; Thompson et al. 2017).

First, the correlation between γ and BIS SPAN is evaluated. In the left panel of Fig. 10, we see that the relationship between γ and the BIS SPAN is linear, and has a slope equal to 0.00072 and a strong Pearson correlation coefficient of R = 0.954. This strong correlation suggests that both γ and BIS SPAN are measuring similar asymmetries for the CCFs. It also provides a conversion between the dimensionless γ parameter into m s⁻¹ using the slope of 0.00072 m s⁻¹.

The right plot of Fig. 10 shows the comparison between the RVs estimated using the SN density and the Normal density. The amplitude of the activity signal is slightly stronger for SN mean RV (in the top right plot the black circles of SN mean RV tend to show more variability), while the signal measured using N mean RV or SN median RV are comparable.

Similar to the analyses presented in Sect. 4, in Fig. 11 we compare the correlations between the asymmetry or the widthparameters of the CCF and the RV. For this analysis, we also include the asymmetry parameters derived in Boisse et al. (2011), V_span and in Figueira et al. (2013), BIS-, BIS+, Bi Gauss, and V_asy, as these authors found those asymmetry parameters more correlated with the RVs than BIS SPAN. It is clear in the case of Alpha Centauri B that the correlation found between γ and SN mean RV is the strongest. The Pearson correlation coefficient is R = 0.74, while the next strongest is R = 0.42 for all the other asymmetry-N mean RV correlations. The correlations between the width and the RV estimates for Alpha Centauri B is also the strongest for the SN parameters, where R = 0.82 for SN FWHM-SN mean RV compared to R = 0.70 for FWHM-N mean RV.

When comparing the correlation between the different Normal and SN parameters in the case of the real data of Alpha Centauri B (see Fig. 11) with the correlations obtained in the SOAP 2.0 simulations (see Sect. 4.3), we observe some significant differences. The correlations between the different parametrizations of the CCF asymmetry and barycenter do not match between the real and simulated data. In the real case, the correlations between γ and SN mean RV, γ and SN median RV and BIS SPAN and N mean RV are all positive. In the case of the SOAP 2.0 simulated data for a facula or a spot, we always find negative correlations. It is therefore not possible to reproduce with SOAP 2.0 the CCF asymmetry variations observed in real observations. On the contrary, the correlations between the CCF width and barycenter match between the real data and a SOAP 2.0 simulated facula. It seems therefore that SOAP 2.0 simulation are able to correctly model the width variation of the CCF, however not its asymmetry. This is probably because in the SOAP 2.0 simulation, a facula is modeled using the same spectrum as a spot with only a different flux. It is well known that facula have a different temperature than spots, and therefore a spectrum that significantly differs (e.g., Cavallini et al. 1985). Looking at other stellar activity simulation like StarSim (Herrero et al. 2016), it seems that simulating a positive correlation between the CCF asymmetry and barycenter is not possible with current tools, and some progress still needs to be made.

Results illustrating the performance of the stellar activity correction proposed in Sect. 3 are indicated in Fig. 12. For Alpha Centauri B, the RV estimated with SN mean RV has a rms that is 35% larger than the rms of the RV estimated with the N mean RV, and the rms of SN median RV is 9% larger than that of the N mean RV. Even though we see these differences in the estimated RV, once we correct for stellar activity using Eq. (5), the rms of the residuals are essentially the same for all three approaches. Although the correlations between the different parameters from the SN density are more sensitive to stellar activity than those obtained with a Normal density fit, the proposed linear model that corrects for stellar activity does not necessarily perform better in the SN case than in the Normal case. The new correction for stellar activity proposed in Sect. 3 performed only slightly better than the usual correction that uses only a linear combination of the width and asymmetry of the CCF.

The results of the statistical tests of the different parameters used for correcting activity can be found in Table 5. The BIS SPAN (coefficient β₂) is not statistically significant for the parameters derived from the Normal density fit. However, all the other parameters in the Normal and SN cases are statistically significant for modeling stellar activity. By analyzing the values of the coefficient of determination, R², we see that the model for SN mean RV is able to capture the highest percentage of variability in the estimated RV. This is not a surprising result since the three different RV estimates have the same RV residual rms after correction for activity, but before correction, SN mean RV had the largest RV rms (see Fig. 12).

When looking at the results discussed in this section, it is likely that the activity signal of Alpha Centauri B is because of faculae. Like observed for the simulated facula in Sect. 4.1, the amplitude of the activity signal is slightly stronger for SN mean RV than for N mean RV orSN median RV; the amplitude of the two latest being comparable. In addition, when applying the proposed correction for activity in the case of the Alpha Centauri B data, the interaction term is significant, which was only the case for the simulated facula in Sect. 4.1. Those are arguments strengthen the findings of Dumusque (2014) who also found evidence for faculae dominating the RV stellar signal of Alpha Centauri B.

Fig. 8 Radial velocity estimates for N mean RV (red dashed line), SN mean RV (black line), or SN median RV (cyan dot-dashed line). In this case, the CCFs have been generated using SOAP 2.0, considering an equatorial 1% spot on the simulated Sun in additionto a planet with a period of one-third of the rotational period of the star and an amplitude of 10 m s−1. The star does one full rotation between phase −0.5 and 0.5; the spot is seen face-on at phase 0.

Open with DEXTER

Fig. 9 Evaluation of the correlation between the RVs and the asymmetry parameters of the simulated data with a 1% spot and an injected planetary signal. The shape of the CCF changes as the spot moves, producing statistically significant correlations only between the estimated RVs and the asymmetry parameter. The correlations between the estimated RVs and the width parameter of the CCF are weaker than the case with only a spot.

Open with DEXTER

5.2 Comparison for HD 192310, HD 10700, HD 215152, and CoRoT-7 of the different CCF parameters derived with the Normal and Skew Normal

In the previous section we evaluated for Alpha Centauri B the improvement obtained by the SN parameters compared to the Normal parameters and the BIS SPAN. We carry out similar analyses for four other main-sequence stars: HD 192310 (K2V; Pepe et al. 2011), HD 10700 (G8V; Feng et al. 2017b), HD 215152 (K3V; Delisle et al. 2018), and finally CoRoT-7 (K0V; Haywood et al. 2014). The same correlation and residual plots shown in the previous section for Alpha Centauri B can be found for those new four stars in Appendix A.

The correlations between the parameters of these additional stars are similar to those obtained for Alpha Centauri B. The correlation between γ and SN mean RV is the strongest among all the asymmetry-RV correlations. Between the width parameters and the estimated RV, the strongest correlation often is between SN FWHM and SN mean RV. However, there is one exception in the case of HD 10700 where the Pearson correlation coefficient between FWHM and N mean RV is equal to R = 0.53, while it is R = 0.42 between SN FWHM and SN mean RV, and R = 0.5 between SN FWHM and SN median RV. As in the case of Alpha Centauri B, positive correlations are always observed between the asymmetry and barycenter of the CCF. This cannot be explained by SOAP 2.0 and to our knowledge by other stellar activity simulators, showing the limit of these tools.

Except for the special case discussed above for HD 10700, the analyses of those four stars, in addition to the analyses on Alpha Centauri B, show that the parameters derived when using a SN density are generally more sensitive to activity. Therefore using the SN parameters, and in particular estimating RV using SN mean RV, can result in better detection of stellar activity than the Normal parameters. More specifically, this is the case for the evaluation of the asymmetry-RV correlations for Alpha Centauri B, HD 10700, HD 215152, HD 192310, and CoRoT-7, and the width-RV correlation for Alpha Centauri B, HD 215152, HD 192310, and CoRoT-7 (see Appendix A).

When correcting for stellar activity for Alpha Centauri B, although the uncorrected RV rms was larger for SN mean RV (compared to the RVs obtained using N mean RV), once corrected for activity using the new model proposed in Sect. 3, both RVs estimates had similar residuals. For HD 10700, HD 215152, and HD 192310, the proposed and usual models were giving similar RV residual rms. However, for CoRoT-7, the new correction is able to provide RV residual rms 23 cm s⁻¹ smaller than that obtained with the usual correction.

Fig. 10 Left panel: correlation between γ and the BISSPAN for Alpha Centauri B. The strong correlation suggests these two parameters are similarly measuring the asymmetry. Top right panel: RVs as functions of Julian Day for Alpha Centauri B in 2010. The RVs are estimated using the mean of a Normal fitted to the CCF (red triangles), or the mean (black circles) or median (cyan pluses) of a SN density fitted to the CCF. Bottom right panel: differences between the RVs estimated with the Normal density and those from the SN density.

Open with DEXTER

5.3 Detection limits when using the estimated RVs from the Normal or Skew Normal models

In the previous section, we saw that the estimated RV resulted in different amplitudes when considering a SN or a Normal density, especially when using SN mean RV. However, once corrected for stellar activity using the linear combination presented in Eq. (5), as shown in the bottom plots of Fig. 12, the rms of the residuals is essentially the same for all three approaches. In this section, we investigate the ability of the three different RV estimators (N mean RV, SN mean RV, and SN median RV) to detect planetary signals among stellar activity. We also compare these RV estimators when using the usual stellar activity correction with the proposed stellar activity model of Eq. (5). To carry out this test, the minimum detected amplitude of an injected planetary signal is estimated at different orbital periods when considering data affected by stellar activity.

In order to obtain CCFs affected by realistic stellar activity signals, the CCFs from Alpha Centauri B used previously were considered. To simulate a planetary signal, the CCFs were blue- or redshifted with the desired amplitude, period, and phase. Several RV data sets with the same stellar signal, but different planetary signals were generated using parameters corresponding to the following grid:

• period of 3, 5, 7, 9, 11, 15, 20, 25, and 30 days,

• amplitude from 0.5 to 3 m s⁻¹ by steps of 0.05 m s⁻¹,

• 10 different phases, evenly sampled between 0 and 2π.

For each of the 4590 simulations we computed the three estimates of RV, namely N mean RV, SN mean RV, and SN median RV. On each set of RV estimates, we performed an analysis similar to Sect. 4.3, i.e., fitting the activity signal using Eq. (5) or the usual correction along with a circular planetary signal (see Eq. (6)). Because of the nonlinearity of the model that includes a planet, a nonlinear least squares algorithm was used for the fit (Levenberg 1944; Marquardt 1963; Teunissen 1990). Such a model requires initial conditions close to the real solution, otherwise the algorithm can converge to a local minimum. Because our goal is to compare the planetary detection limits using the three different RV estimates and the two different activity models proposed, and not to discuss what is the best method to explore the parameter space, we initialized the minimization algorithm to the real period of the planetary signal injected to avoid getting stuck in a local minimum. We also selected as initial amplitude the peak-to-peak amplitude of the estimated RVs. The argument of periapsis t₀ was initialized to the time when the RV was crossing 0 since we use a sinusoidal function to fit the planetary signal (see Eq. (7)).

Once the parameters involved in Eq. (6) were estimated, signals in the residuals, defined as RV −RV_activity, were analyzed using a generalized Lomb–Scargle periodogram (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009). If a signal with a P-value⁵ smaller than 1% had a period compatible with the injected planetary period within an error budget of 20%, the signal was considered significant and the corresponding planet considered detected. For each period considered, we searched for the minimum amplitude at which at least 80% of the planets with different phases were detected. This minimum amplitude detected as a function of period is shown in Fig. 13 for the three different RV estimates (N mean RV, SN mean RV, and SN median RV) when using the new stellar activity correction proposed in this paper (see Eq. (5)), and when using the usual activity correction. We can see that the new correction for stellar activity based on Eq. (5) improves the detection limit of the exoplanet by 12% on average compared to the usual approach, and the three estimators of RV give similar detection limits. These results therefore suggest that any of the RV estimators can be used when searching for a planetary signal in RV data contaminated by stellar activity, and using our new model to account for stellar activity allows us to detect planetary signals with a slightly smaller amplitude than the usual correction that uses only a linear correlation with the FWHM and BIS SPAN.

	Fig. 11 Topthree rows: correlations between the asymmetry parameters and their corresponding estimated RVs for Alpha Centauri B. Bottom row: correlation between the FWHM and the estimated RVs. The correlationsare stronger when using parameters derived from the SN fit than the Normal fit. The estimated R’s are all statistically significant.
Open with DEXTER

Fig. 12 Top panels: RVs (black dots) for Alpha Centauri B estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (new corr. std, black dots) and the residuals from the usual correction (usual corr. std, blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals have a smaller systematic component when using the proposed model of Eq. (5) (black dots) compared to the usual model (blue triangles).

Open with DEXTER

6 Estimation of standard errors for the CCF parameters

In this section, we investigate how the photon noise influences the CCF parameters derived either by a Normal density or SN density fit. Because a CCF is obtained from a cross-correlation, each point of a CCF is correlated with the other points. Therefore, we cannot simply vary each point in the CCF by their respective error bars and then recalculate the best SN or Normal density fit to see how the CCF noise influences the estimation of the parameters of interest (i.e., N mean RV, SN mean RV, SN median RV, FWHM, SN FWHM, BIS SPAN, and γ). Instead, we focus on the individual spectrum where each individual point can be considered independent from the others. The standard error on each point of a spectrum is given by the photon noise, which follows a Poisson distribution and is therefore estimated by taking the square root of the measured flux.

The following method was carried out in order to estimate the error bars on the different parameters derived from the CCF. We first modify the values of all the points in the spectrum given their respective error bars. To do so random Gaussian noise with standard deviation the square root of the flux was added across each spectrum. The CCF was calculated using this spectrum according to the method presented in Pepe et al. (2002), then fit by either a Normal or SN density with the parameters recorded. This process was repeated a hundred times to obtain a distribution for each CCF parameter, and the standard deviations of the resulting distributions provide estimates of the standard errors for the CCF parameters.

The standard errors were computed for each CCF parameter for the HARPS measurements of HD 215152, HD 192310, and CoRoT-7. These three stars include measurements that cover the range of S/N measured at 550 nm (S/N550) from 10 to 500, which represent the very low S/N limit and saturation limit of the HARPS detector, respectively. HD 10700 and Alpha Centauri B were not included because they have a large number of measurements, which would require a substantial computational effort. The variation of the noise for each CCF parameter as a function of S/N550 is shown in Fig. 14. The top row shows the standard errors of the three different estimated RVs, the width, and the asymmetry estimates. Because BIS SPAN and γ do not have the same units, the estimated slopes of the correlation between those two parameters to transform γ in m s⁻¹ were used (see Fig. 10 and Table A.1 for the value of the slope for each star). The bottom row shows the ratio between the standard errors measured when using the SN parameters and the Normal parameters. Values smaller (larger) than one imply that standard errors from the SN parameters are more (less) precise than the Normal parameters.

The standard errors for the different RV estimates all appear to follow a similar exponential decay as a function of S/N, even though the measurements are from three different stars. This suggests that, for the considered stars, the precision in RV is mainly driven by the S/N of the analyzed spectra. As shown in Bouchy et al. (2005), the RV precision is proportional to the S/N, the FWHM of the CCF, and its contrast. In our case all three studied stars are main-sequence K dwarfs, which implies that their CCF FWHM and contrast are similar and explains why the RV precision is driven by the S/N only.

When comparing the three different estimates for the RV, SN mean RV has standard errors that are 60% larger than N mean RV. However, SN median RV gives errors 10% more precise than N mean RV. The parameters describing the width of the CCF, FWHM, and SN FWHM have comparable standard errors. Finally, for the asymmetry parameters, γ has standard errors that are 15% more precise than BIS SPAN. In conclusion, when fitting a SN density to the CCF and when using SN median RV as the RV estimate, we are able to improve the precision on the estimated RV by 10%. Using the SN density, we are also able to improve by 15% the precision on the estimated asymmetry parameter of the CCF. However, SN mean RV should not be used to derive precise RV estimates as the precision on this parameter is 60% worse than the precision on the RVs derived from N mean RV.

Fig. 13 Detection limits of planetary signals once the stellar activity signal is removed from the raw RVs using the model proposed in Eq. (5) (solid lines) or the usual correction based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit (dashed lines). The correction for stellar activity based on Eq. (5) improves on average the detection limit by 12% and the different RV estimators have similar detection limits.

Open with DEXTER

7 Discussion

When fitting a SN density shape to the CCF, parameters used to estimate the RV, defined as the CCF barycenter, the amplitude, sometimes called the CCF contrast, the width and asymmetry of the CCF can all be estimated in a single model framework. For the estimation of the RV, we investigated the use of the mean and median of the SN density. The width is derived using the variance of the SN density (SN ) and the asymmetry by using γ, Skewness parameter of the SN density.

To evaluate the performance of the proposed SN framework, tests on both simulated and real data were carried out and compared to the commonly employed approach of fitting a Normal density shape to the CCF to get access to the RV and FWHM, and then separately deriving the BIS SPAN to estimate the CCF asymmetry. The simulated CCFs were generated using the SOAP 2.0 code, which can simulate activity signals induced by a spot or a facula on a solar-like star. To simulate realistic data, we considered a S/N of 100, which is typical for high-precision RV observations.

The results of the simulation study suggest that at least one of the parameters derived from the SN density fit is equally or more sensitive to activity than the parameters obtained by the usual Normal method, making this or these parameters more useful indicators of activity. Sensitivity was measured by the strength of the correlation between thedifferent SN or Normal parameters. In the case of a spot, the strongest correlations are found between γ and SN median RV and betweenBIS SPAN and N mean RV, therefore making the SN parameters equally sensitive to activity. For the facula case, the strongest correlation is between SN FWHM and SN mean RV with a correlation coefficient of R = 0.92. The correlation between the parameters derived from the Normal fit are much weaker with correlation coefficients between BIS SPAN and N mean RV of R = −0.61 and betweenFWHM and N mean RV of R = 0.47. The SN parameters continued to have stronger correlations than the Normal parameters in the setting where a planetary signal was added the SOAP 2.0 with a single spot.

Looking at real data, we arrive at a similar conclusion that the SN parameters are more sensitive to activity. While real data confirms that the correlation between CCF width and barycenter is always stronger for the SN parameters, they also show that it is the case for the correlations between the CCF asymmetry and barycenter. In this later case, all stars studied in this paper show a positive correlation, which cannot be explained by SOAP 2.0, as the spot and facula simulations only show negative correlations between the CCF asymmetry and barycenter. This discrepancy between simulations and real data could be because SOAP 2.0 uses the spectrum of a spot as input to model the activity induced by a facula. Because the temperature between a spot and a facula is significantly different, their spectra should be different. Additionally, there are expected to be multiple active regions on a star at different locations in longitude and latitude, while the SOAP 2.0 data used in the simulation study included only a single active region on the equator in order to isolate the effects of those active regions.

In the real cases, the parameters derived from the SN are always more sensitive to activity than the parameters derived from the Normal. There is only one exception in the case of the width-barycenter correlation for HD 10700, however, the difference in correlation between the SN and Normal parameters is rather small, from R = 0.42 to 0.53, respectively. Also, the correlation between the asymmetry and SN mean RV is consistently stronger than the parametrization of the CCF presented in Boisse et al. (2011) and Figueira et al. (2013). Because the apparent RV signal induced by activity results in a stronger correlation with the SN parameters than between the apparent RV signal and the FWHM of the CCF or its BIS SPAN, this suggests the SN model of the CCF can lead to a better understanding of the spurious variations in RV caused by stellar activity.

Considering the different RV estimates of the real data, the amplitude of stellar activity tends to be largest for SN mean RV, followed by SN median RV, and N mean RV, which behave similarly. This implies that the mean of the SN density appears to be more sensitive to variation in the CCF shape than the median of the SN or the mean of the Normal.

Having an estimator of RV that is more sensitive to stellar activity, such as SN mean RV, can also help to better probe stellar rotational periods or to better understand the covariance of stellar signals when fitting a Gaussian process to the RVs (e.g., Faria et al. 2016; Haywood et al. 2014). We saw that the SN mean RV estimator is 60% noisier than the N mean RV estimator. However, when looking at bright stars such as α Centauri B or HD 10700, increasing the photon noise by 60% does not have a significant impact on the RV precision because the instrumental noise dominates the data. Therefore, for bright targets, stellar activity can be better characterized by using SN mean RV as this RV estimate is more sensitive to it.

We also propose a new model to correct the estimated RV data for stellar activity signals. Generally, when fitting for planetary signals, it is common to use a model composed of one or several Keplerian signals to account for the planets, in addition to a linear combination of the FWHM and BIS SPAN to account for stellar activity signals. The proposed model adds a term to the linear model to account for the amplitude of the CCF and an interaction term between the estimated asymmetry and the width parameters. Using the simulated data from SOAP 2.0, this new model reduces the effect of the stellar activity signal by factors of about 2 and 3.5 over the usual model for the facula and spot, respectively.

Even if the different RV estimators derived by the Normal or SN fit result in different amplitudes, once the proposed correction for stellar activity is applied the residuals of the model have similar rms. When comparing the activity correction proposed in this paper with the usual correction that only uses a linear combination of the CCF asymmetry and width, the new proposed correction almost entirely explains the spurious variations in RV for the simulations based on the presence of a facula or a spot. However, when moving to real data, there is just a slight improvement by using the proposed correction function for stellar activity. Additional analysis can be performed for new data sets to see if certain components of the model proposed in this work are not relevant and, therefore, could be removed.

A test was carried out to see if some RV estimates were better at finding planets in RV data affected by observed stellar signals. The new correction based on Eq. (5) proposed in this paper to mitigate the effect from stellar activity slightly improves the detection limit with respect to the usual limit based on RV_activity = β₀ + β₁γ + β₂SN FWHM for the SN fit and on RV_activity = β₀ + β₁BIS SPAN + β₂FWHM for the Normal fit. Concerning the definition of the RV using the SN or the Normal fit, all three of the different RV estimators give similar detection limits. Therefore it seems that any of the RV estimators can be used to search for planetary signals.

Finally, we investigated the precision of each of the SN and Normal parameters including the BIS SPAN. It turns out that SN mean RV should not be used to get precise RVs as the standard errors on this parameter is 60% greater than for N mean RV. However, SN median RV is 10% more precisethan N mean RV. Regarding the asymmetry estimates, we observe that γ has a precision 15% better than the BIS SPAN.

Fig. 14 Results of the bootstrap analyses on the stars HD 215152, HD 192310, and CoRoT-7. Top panels: comparison between the standard errors from the bootstrap analysis of the estimated RVs, FWHM, and asymmetry parameters using the SN fit and the common strategy (Normal fit and BIS SPAN). Bottom panels: ratio between the standard errors estimated on the parameters derived from the common strategy and the corresponding standard errors estimated on the parameters derived from the SN fit. When using SN mean RV (black circles), the standard errors are on average 60% larger than the standard errors of N mean RV (red triangles). However, the standard errors for SN median RV (cyan pluses) are on average 10% smaller than the standard errors coming from the N mean RV. The use of the asymmetry SN parameter γ leads to standard errors on average 15% smaller than the standard errors related to the BIS SPAN. We note that for asymmetry, the error in BIS SPAN is in m s−1. To be able tocompare the errors in γ and BIS SPAN, we multiplied the error in γ by the slope of the correlation between γ and BIS SPAN.

Open with DEXTER

8 Conclusions

When searching for low-mass exoplanets using the RV technique, it is necessary to retrieve precise estimates of the RV and also to account for variations induced by stellar activity in order to avoid false detections. Stellar activity such as spots and faculae can lead to shape variations in the spectra features, which then results in shape variations of the CCF. The correlations between the width or asymmetry of the CCF and the estimated RV are commonly used as a way to detect if the RVs are affected by stellar activity signals. Because the presence of real planets would result in only a shift in the CCF (not a change of its shape), strong correlations between the shape features of the CCF and the estimated RVs suggest that stellar activity may be present.

In this paper, a new approach for quantifying shape changes in the CCF is proposed using the SN density, which can be used to estimate with a single fit the RV, width, and Skewness of the CCF. This new method is compared to a commonly used method based on a Normal density fit to the CCF. The mean of the Normal density is used as the estimated RV and the FWHM estimates the width of the CCF. Because the Normal density does not have any Skewness, another method is necessary to estimate the asymmetry of the CCF, such as the often employed BIS SPAN. In addition, the proposed SN approach is compared to other parameterizations of the CCF asymmetry, which have been shown to be sensitive to activity signals (Boisse et al. 2011; Figueira et al. 2013).

In the different tests carried out for this work, the SN parameters performed at least as well as, and most of the times better than, the parameters from the Normal approach and BIS SPAN. The SN parameters γ, SN FWHM, and SN mean RV consistently had stronger correlations than those between any of the parameters derived by the Normal and BIS SPAN, or the different asymmetry parametrizations presented in Boisse et al. (2011) and Figueira et al. (2013). This suggests the SN parameters may be better at probing stellar activity signals than the other methods. In addition, the uncertainties measured on SN median RV and γ are 10% and 15%, respectively, smaller than the corresponding uncertainties on N mean RV and BIS SPAN, although SN means RV had uncertainties 60% greater than N mean RV.

Because of the advantages of using the proposed SN approach over the commonly employed approach based on the Normal density fit to the CCF and the BIS SPAN or the asymmetry parameters described in Boisse et al. (2011) and Figueira et al. (2013), a SN density model for the CCF may be more useful for detecting stellar activity than the previously proposed parametrizations. Correlations between γ and SN mean RV and between the width and SN mean RV can be used to probe stellar activity signals in RV data, and SN median RV can be used to estimate RV. We also proposed a new model to correct the estimated RV data for stellar activity signals, by using the amplitude of the CCF and an interaction term between the estimated asymmetry and the width parameters. Using simulated data from SOAP 2.0, this new proposed correction reduces the effect of the stellar activity signal by an additional 14.5–15% and 5.3–5.8% over the usual model for facula and spot, respectively. When applying this model on real data, we observe that planetary detection limits are improved by a non-negligible 12%.

Acknowledgements

The authors thank Yale’s Center for Research Computing for their help and resources with some of the computational aspects of this work. X.D. is grateful to The Branco Weiss Fellowship–Society in Science for its financial support. J.C.K. was partially supported by the National Science Foundation under Grant AST 1616086 and by the National Aeronautics and Space Administration under grant 80NSSC18K0443. U.S. was partially supported by Fondazione CARIPARO and thanks the IT-University of Helsinki for the computational resources provided to execute part of the analyses of the present work. The authors are grateful to all technical and scientific collaborators of the HARPS Consortium, ESO Headquarters, and ESO La Silla who have contributed with their extraordinary passion and valuable work to the success of the HARPS project.

Appendix A Additional table and figures

In this appendix, a similar analysis as that presented in Sect. 5 is discussed for four main-sequence stars: HD 192310 (K2V; Pepe et al. 2011), HD 10700 (G8V; Feng et al. 2017b), HD 215152 (K3V; Delisle et al. 2018), and finally CoRoT-7 (K0V; Haywood et al. 2014). The latest HARPS data for these stars can be found on the ESO archive.

Table A.1 summarizes the results obtained by the SN and Normal density models of the CCF. Theseresults are consistent with those from the analysis of Alpha Centauri B. The correlation between γ and SN mean RV is stronger than the correlation between the BIS SPAN and N mean RV orbetween the asymmetry parameters derived in Boisse et al. (2009) and Figueira et al. (2013) and N mean RV for all the considered stars. The correlation between SN FWHM and SN mean RV is stronger than the correlation between FWHM and

N mean RV for three of the four stars. Also for all these stars, the originally estimated RVs were corrected from spurious variations caused by stellar activity using Eq. (5), and Figs. A.2, A.4, A.6, and A.8 show the corrected RVs. Once corrected from stellar activity, the Normal and SN residuals are comparable for the stars HD 192310, HD 10700, and HD 215152. However, the rms of the residuals for CoRoT-7are 0.23 m s⁻¹ lower for the SN model than the Normal model. The average S/N at 550 nm for the stars HD 10700, HD 192310, HD 215152, and CoRoT-7 are, respectively, 273, 207, 141, and 44. CoRoT-7 has therefore on average a much lower S/N at 550 nm than the others stars, which could be a potential explanation for this small improvement. Additional tests should be performed to confirm this statement.

	Fig. A.1 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for HD 192310. Bottom row: correlations between the FWHM and the estimated RVs. The correlations are consistently stronger when using parameters derived from the SN than the Normal. The estimated R are all statistically significant.
Open with DEXTER

Fig. A.2 Toppanels: RVs (black dots) for HD 192310 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.

Open with DEXTER

Fig. A.3 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for HD 10700. Bottom row:correlations between the FWHM and RVs for HD 10700. The correlations are consistently stronger when using SN mean RV compared to N mean RV for the asymmetry parameters; however, the correlation between the FWHM and the N mean RV, only for this quiet star, is stronger the analogous correlations with the estimated SN RVs. The estimated R are statistically significant, except for the correlation between FIG BIS and RV (p-values = 0.36).

Open with DEXTER

Fig. A.4 Toppanels: RVs (black dots) for HD 10700 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.

Open with DEXTER

Fig. A.5 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for HD 215152. Bottom row:correlations between the FWHM and the RVs for HD 215152. The correlations are consistently stronger when using SN mean RV compared to N mean RV. The p-values associated with each R are not statistically significant for the correlation between N mean RV and BIS SPAN (p-values = 0.27), the correlation between N mean RV and FIG BIS- (p-values = 0.05), the correlation between SN median RV and SN FWHM (p-values = 0.5), and the correlation between N mean RV and FWHM (p-values = 0.2).

Open with DEXTER

Fig. A.6 Toppanels: RVs (black dots) for HD 215152 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.

Open with DEXTER

Fig. A.7 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for CoRoT7. Bottom row:correlations between the FWHM and the RVs for CoRoT 7. The correlations are consistently stronger when using parameters derived from the SN than the Normal. The p-values associated with each R are not statistically significant for the correlation between N mean RV and BIS SPAN (p-values = 0.23) and the correlation between N mean RV and FIG BIS- (p-values = 0.11).

Open with DEXTER

Fig. A.8 Toppanels: RVs (black dots) for CoRoT 7 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals have a smaller systematic component when using the proposed model of Eq. (5) (black dots) compared to the usual model (blue triangles). Moreover, once corrected for stellar activity using Eq. (5), the remaining standard deviation from the SN models are 0.334 m s−1 smaller than the remaining standard deviation of the Normal model.

Open with DEXTER

References

All Tables

All Figures

	Fig. 1 Density function of a random variable Y following a SN distribution SN(ξ, ω, α) with location parameter ξ = 0, scale parameter ω = 1 and differentvalues of the Skewness parameter α indicated by different colors and line types. We note that the solid black line has an α = 0 making it a Normal distribution.
Open with DEXTER
In the text

	Fig. 2 RV estimates for N mean RV (red dashed line), SN mean RV (black line), and SN median RV (cyan dotted-dashed line). In this case, the CCFs were generated using SOAP 2.0 with an equatorial3% facula on the simulated Sun. The star does one full rotation between phase −0.5 and 0.5; the facula are seen face-on for phase 0. The variations observed in SN mean RV are notably different from the variations measured in SN median RV and N mean RV.
Open with DEXTER
In the text

	Fig. 3 Left panels: correlations between the different asymmetry parameters and their corresponding RV estimates in the case of an equatorial 3% facula on the simulated Sun. Right panels: correlations between the different width parameters and their corresponding RV estimates for the same facula. In the presence of a facula, both the shape and the width of the CCF change as the star rotates, producing statistically significant correlations for all the cases except for the correlation between SN mean RV and SN GAMMA (P-value = 0.27).
Open with DEXTER
In the text

	Fig. 4 Top panels: spurious estimated RVs (black dots) caused by a facula in the simulated data using a Normal and a SN model, the estimated RVs using Eq. (5) (red pluses), and the estimated RVs using the usual correction for stellar activity (green triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (red pluses) and the residuals from the usual correction (green triangles). The standard deviations are also reported in the legend, and the residuals have a smaller systematic component when using the proposed model of Eq. (5) compared to the usual model. The tests of statistical significance on the parameters are presented in Table 2.
Open with DEXTER
In the text

	Fig. 5 RV estimates for N mean RV (red dashed line), SN mean RV (black line) or SN median RV (cyan dotted-dashed line) using CCFs generated from SOAP 2.0 with an equatorial 1% spot on the simulated Sun. The star does one full rotation between phase −0.5 and 0.5; the spot is seen face-on at phase 0. The SN mean RV seems to have the smallest spurious variations caused by the spot.
Open with DEXTER
In the text

	Fig. 6 Left panels: correlations between the different asymmetry parameters and their corresponding RV estimates in the case of an equatorial 1% spot on the simulated Sun. Right panels: correlations between the different width parameters and their corresponding RV estimates for the same spot. In the presence of a spot, both the shape and width of the CCF change as the star rotates. However, only the asymmetry produces a statistically significant correlation with the different RV estimates. The width parameters and their corresponding RV estimates present weak correlations and, in general, much weaker correlations compared to the results obtained when an equatorial 3% facula is present on the simulated Sun.
Open with DEXTER
In the text

	Fig. 7 Top panels: spurious estimated RVs (black dots) caused by a spot in the simulated data, the estimated RVs using Eq. (5) (red pluses), and the estimated RVs using the usual correction for stellar activity (green triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (red pluses) and the residuals from the usual correction (green triangles). The standard deviations are also reported in the legend, and the residuals have a smaller systematic component when using the proposed model compared to the usual model. The tests of statisticalsignificance on the parameters are presented in Table 3.
Open with DEXTER
In the text

	Fig. 8 Radial velocity estimates for N mean RV (red dashed line), SN mean RV (black line), or SN median RV (cyan dot-dashed line). In this case, the CCFs have been generated using SOAP 2.0, considering an equatorial 1% spot on the simulated Sun in additionto a planet with a period of one-third of the rotational period of the star and an amplitude of 10 m s−1. The star does one full rotation between phase −0.5 and 0.5; the spot is seen face-on at phase 0.
Open with DEXTER
In the text

	Fig. 9 Evaluation of the correlation between the RVs and the asymmetry parameters of the simulated data with a 1% spot and an injected planetary signal. The shape of the CCF changes as the spot moves, producing statistically significant correlations only between the estimated RVs and the asymmetry parameter. The correlations between the estimated RVs and the width parameter of the CCF are weaker than the case with only a spot.
Open with DEXTER
In the text

	Fig. 10 Left panel: correlation between γ and the BISSPAN for Alpha Centauri B. The strong correlation suggests these two parameters are similarly measuring the asymmetry. Top right panel: RVs as functions of Julian Day for Alpha Centauri B in 2010. The RVs are estimated using the mean of a Normal fitted to the CCF (red triangles), or the mean (black circles) or median (cyan pluses) of a SN density fitted to the CCF. Bottom right panel: differences between the RVs estimated with the Normal density and those from the SN density.
Open with DEXTER
In the text

	Fig. 11 Topthree rows: correlations between the asymmetry parameters and their corresponding estimated RVs for Alpha Centauri B. Bottom row: correlation between the FWHM and the estimated RVs. The correlationsare stronger when using parameters derived from the SN fit than the Normal fit. The estimated R’s are all statistically significant.
Open with DEXTER
In the text

	Fig. 12 Top panels: RVs (black dots) for Alpha Centauri B estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (new corr. std, black dots) and the residuals from the usual correction (usual corr. std, blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals have a smaller systematic component when using the proposed model of Eq. (5) (black dots) compared to the usual model (blue triangles).
Open with DEXTER
In the text

	Fig. 13 Detection limits of planetary signals once the stellar activity signal is removed from the raw RVs using the model proposed in Eq. (5) (solid lines) or the usual correction based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit (dashed lines). The correction for stellar activity based on Eq. (5) improves on average the detection limit by 12% and the different RV estimators have similar detection limits.
Open with DEXTER
In the text

	Fig. 14 Results of the bootstrap analyses on the stars HD 215152, HD 192310, and CoRoT-7. Top panels: comparison between the standard errors from the bootstrap analysis of the estimated RVs, FWHM, and asymmetry parameters using the SN fit and the common strategy (Normal fit and BIS SPAN). Bottom panels: ratio between the standard errors estimated on the parameters derived from the common strategy and the corresponding standard errors estimated on the parameters derived from the SN fit. When using SN mean RV (black circles), the standard errors are on average 60% larger than the standard errors of N mean RV (red triangles). However, the standard errors for SN median RV (cyan pluses) are on average 10% smaller than the standard errors coming from the N mean RV. The use of the asymmetry SN parameter γ leads to standard errors on average 15% smaller than the standard errors related to the BIS SPAN. We note that for asymmetry, the error in BIS SPAN is in m s−1. To be able tocompare the errors in γ and BIS SPAN, we multiplied the error in γ by the slope of the correlation between γ and BIS SPAN.
Open with DEXTER
In the text

	Fig. A.1 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for HD 192310. Bottom row: correlations between the FWHM and the estimated RVs. The correlations are consistently stronger when using parameters derived from the SN than the Normal. The estimated R are all statistically significant.
Open with DEXTER
In the text

	Fig. A.2 Toppanels: RVs (black dots) for HD 192310 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.
Open with DEXTER
In the text

	Fig. A.3 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for HD 10700. Bottom row:correlations between the FWHM and RVs for HD 10700. The correlations are consistently stronger when using SN mean RV compared to N mean RV for the asymmetry parameters; however, the correlation between the FWHM and the N mean RV, only for this quiet star, is stronger the analogous correlations with the estimated SN RVs. The estimated R are statistically significant, except for the correlation between FIG BIS and RV (p-values = 0.36).
Open with DEXTER
In the text

	Fig. A.4 Toppanels: RVs (black dots) for HD 10700 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.
Open with DEXTER
In the text

	Fig. A.5 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for HD 215152. Bottom row:correlations between the FWHM and the RVs for HD 215152. The correlations are consistently stronger when using SN mean RV compared to N mean RV. The p-values associated with each R are not statistically significant for the correlation between N mean RV and BIS SPAN (p-values = 0.27), the correlation between N mean RV and FIG BIS- (p-values = 0.05), the correlation between SN median RV and SN FWHM (p-values = 0.5), and the correlation between N mean RV and FWHM (p-values = 0.2).
Open with DEXTER
In the text

	Fig. A.6 Toppanels: RVs (black dots) for HD 215152 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.
Open with DEXTER
In the text

	Fig. A.7 Top three rows: correlations between the asymmetry parameters and their corresponding RVs for CoRoT7. Bottom row:correlations between the FWHM and the RVs for CoRoT 7. The correlations are consistently stronger when using parameters derived from the SN than the Normal. The p-values associated with each R are not statistically significant for the correlation between N mean RV and BIS SPAN (p-values = 0.23) and the correlation between N mean RV and FIG BIS- (p-values = 0.11).
Open with DEXTER
In the text

	Fig. A.8 Toppanels: RVs (black dots) for CoRoT 7 estimated using a SN and a Normal fit. Bottom panels: residuals from the model fit using Eq. (5) (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on RVactivity = β0 + β1γ + β2SN FWHM for the SN fitand on RVactivity = β0 + β1BIS SPAN + β2FWHM for the Normal fit. The residuals have a smaller systematic component when using the proposed model of Eq. (5) (black dots) compared to the usual model (blue triangles). Moreover, once corrected for stellar activity using Eq. (5), the remaining standard deviation from the SN models are 0.334 m s−1 smaller than the remaining standard deviation of the Normal model.
Open with DEXTER
In the text