Measuring precise radial velocities and cross-correlation function line-profile variations using a Skew Normal density

11/30/2018
by   Umberto Simola, et al.
Helsingin yliopisto
0

Stellar activity is one of the primary limitations to the detection of low-mass exoplanets using the radial-velocity (RV) technique. We propose to estimate the variations in shape of the CCF by fitting a Skew Normal (SN) density which, unlike the commonly employed Normal density, includes a skewness parameter to capture the asymmetry of the CCF induced by stellar activity and the convective blueshift. The performances of the proposed method are compared to the commonly employed Normal density using both simulations and real observations, with different levels of activity and signal-to-noise ratio. When considering real observations, the correlation between the RV and the asymmetry of the CCF and between the RV and the width of the CCF are stronger when using the parameters estimated with the SN density rather than the ones obtained with the commonly employed Normal density. Using the proposed SN approach, the uncertainties estimated on the RV defined as the median of the SN are on average 10 Normal. The uncertainties estimated on the asymmetry parameter of the SN are on average 15 Slope Span (BIS SPAN), which is the commonly used parameter to evaluate the asymmetry of the CCF. We also propose a new model to account for stellar activity when fitting a planetary signal to RV data. Based on simple simulations, we were able to demonstrate that this new model improves the planetary detection limits by 12 account for stellar activity. The SN density is a better model than the Normal density for characterizing the CCF since the correlations used to probe stellar activity are stronger and the uncertainties of the RV estimate and the asymmetry of the CCF are both smaller.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 13

page 14

page 16

page 19

page 20

page 21

page 22

11/03/2017

Improving Exoplanet Detection Power: Multivariate Gaussian Process Models for Stellar Activity

The radial velocity method is one of the most successful techniques for ...
03/25/2021

Parameter estimation and model selection for water sorption in a wood fibre material

The sorption curve is an essential feature for the modelling of heat and...
01/26/2020

On truncated multivariate normal priors in constrained parameter spaces

We show that any lower-dimensional marginal density obtained from trunca...
03/30/2020

Thermophysical modelling and parameter estimation of small solar system bodies via data assimilation

Deriving thermophysical properties such as thermal inertia from thermal ...
02/18/2020

Constraining the recent star formation history of galaxies : an Approximate Bayesian Computation approach

[Abridged] Although galaxies are found to follow a tight relation betwee...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

When working with radial-velocities 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. 2017; 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 them, 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 to 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 the 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 the CCF 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 will be observed between the RV and chromospheric activity indicators like (R) or H- (Boisse et al. 2009; Dumusque et al. 2012; Robertson et al. 2014), but also between the RV and the FWHM of the CCF or its BIS SPAN (Queloz et al. 2001; Boisse et al. 2009; Queloz 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 (R), the FWHM, and the BIS SPAN (Dumusque et al. 2017; Feng et al. 2017)

. 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 , which is the difference between the RV measured by fitting a Normal density to the upper and the lower parts of the CCF. This CCF asymmetry parameter is shown to be more sensitive than the BIS SPAN at low signal-to-noise ratio (S/N). Figueira et al. (2013) studied the use of new indicators, BIS-, BIS+, bi-Gauss and . The authors were 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 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 bi-Gauss, the RV and the 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, the FWHM and the asymmetry of the CCF, as this function includes a skewness parameter (Azzalini 1985).

The paper is organized as follows. In Sec. 2 we introduce the SN density, describe its applicability for modelling the CCF and study how the SN parameters relate to the mean of the Normal density, the FWHM, and the BIS SPAN of the CCF. In Sec. 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. 2017). In Sec. 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 Sec. 5. In Sec. 6 error bars are computed for the different estimated CCF parameters. Finally the discussion of the results and the conclusions are presented, respectively, in Secs. 7 and  8.

2 The Skew Normal distribution

The Skew Normal (SN) distribution is a class of probability distributions which 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 Normal distribution’s mean and standard deviation, a third parameter which 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 can be written as

(1)

where and are, respectively, the density function and the distribution function of a standard Normal distribution111A standard Normal distribution is a Normal distribution with a mean of 0 and a standard deviation of 1.. 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 ( and ) are displayed in Fig. 1. A usual Normal distribution is the special case of the SN distribution when the skewness parameter is equal to zero222This can be seen from Eq. def:snd_gen. If then and therefore which is the density of a Normal distribution. Note that because is the the probability that a standard Normal random variable is less than or equal than 0..

Figure 1: Density function of a random variable following a SN distribution with location parameter , scale parameter and different values of the skewness parameter indicated by different colors and line types. Note that the solid black line has an making it a Normal distribution.

For reasons related to the interpretation of the parameters in Eq. def:snd_gen 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. eq:snd_cp, 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 displayed 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 displayed in Fig. 1.

Median
-3 -0.757 0.427 -0.667 -0.672
0 0.000 1.000 0.000 0.000
2 0.714 0.491 0.454 0.655
6 0.787 0.381 0.891 0.674
10 0.794 0.370 0.956 0.674
Table 1: CP values along with the median corresponding to the values shown in Fig. 1, with location parameter and scale parameter . Values are rounded to three decimal places.

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

2.1 Fitting the Skew Normal density to the CCF

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

(3)

where C is 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, respectively, the mean, the variance, and the skewness of the SN as defined above. The values 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 FWHM333FWHM with standard deviation 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 Sec. 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 Sec. 5 and Sec. 6. The width of the SN, SN FWHM, is defined in the same way as for the Normal density444Note that SN FWHM does not correspond to the width of the SN density at half maximum like in the Normal case., 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, , which in its general form is defined as:

(4)

where and are two quantitative variables, indicates the covariance between and , and and represent their standard deviations. A

-value for the statistical test with null hypothesis

is also generally provided.

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 RV signals 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, the 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 the FWHM to take into account the signal induced by stellar activity (e.g. Dumusque et al. 2017; Feng et al. 2017). We propose to add additional parameters in the model to correct for stellar activity: (i) an amplitude parameter A of the CCF (referred to as the CCF contrast) and (ii) 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 and covariance matrix equal to (

defined as the identity matrix). The contrast parameter

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) will be better understood through the results of the examples presented in Sec. 4. This interaction term can account for possible interactions between SN FWHM (or FWHM) and (or BIS SPAN), meaning that each variables’ association with the response, , depends also on the other variable.

The proposed model is analyzed using statistical tests on the parameters , , , and where the null hypothesis is , for . The significance level for the tests are set at . The coefficient of determination, , is used to assess how well the proposed linear combination of variables accounts for the variability of .

The proposed function defined in Eq. eq:RV:correction 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 between the asymmetry and width parameters and therefore the problem of multicollinearity appeared to be mitigated. We also investigated the statistical significant of the interactions between and the width of the CCF, and 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 for modelling the CCF and the benefit of using the proposed correction for stellar activity (see Eq. eq:RV:correction), we begin by considering a simulation study using spectra generated from the Spot Oscillation And Planet 2.0 code (SOAP 2.0, 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 the 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 modelled with a solar disc 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 K and the temperature difference of the facula depends on the centre-to-limb angle , 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. 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.

Figure 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 equatorial 3% facula on the simulated Sun. The star does one full rotation between phase -0.5 and 0.5, with the facula being seen face-on for phase . The variations observed in SN mean RV are notably different from the variations measured in SN median RV and N mean RV.

Correlations between the different RV estimates and the different CCF asymmetry or width estimates are displayed 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, with Pearson correlation coefficient values of =-0.11, -0.55 and -0.61, respectively. Therefore, it seems 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 (). Then follows the correlation between SN FWHM and SN median RV (), and finally FWHM and N mean RV (). Over all the different correlations in the case of a facule, the strongest one 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.

Figure 3: (left) Correlations between the different asymmetry parameters and their corresponding RV estimates in the case of an equatorial 3% facula on the simulated Sun. (right) 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).

The RV variations displayed in Fig. 2 are caused only by stellar activity, in this case a facula. We applied the activity correction proposed in Eq. eq:RV:correction to evaluate the performance of the model in this setting. The results of this correction are displayed in Fig. 4 and the statistical tests on the coefficients involved in Eq. eq:RV:correction 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 larger than . In addition, the proposed linear model allows to reduce the activity effect from a RV rms of 8.02  m sto 1.02  m swhen considering the CCF contrast and FWHM, and the BIS SPAN. When using the parameters derived from the SN, the improvements for SN median RV and SN mean RV are 7.07  m sdown to 0.88  m sand 5.9  m sdown to 0.88  m s, respectively. When comparing the correction proposed in Eq. eq:RV:correction 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.

Figure 4: (top) The 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. eq:RV:correction (red pluses), and the estimated RVs using the usual correction for stellar activity (green triangles), based on for the SN fit and on for the Normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (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. eq:RV:correction compared to the usual model. The tests of statistical significance on the parameters are presented in Table 2.
Parameter N mean RV SN mean RV SN median RV
Table 2: P-values for the estimated coefficients from the model in Eq. eq:RV:correction for correcting stellar activity induced by an equatorial 3% facula on the simulated Sun. 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. The estimated show that the proposed correction for stellar activity explains the vast majority of the spurious variability present in the different RV estimates.

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. Note that this is a large spot for the Sun, as large spots are generally on the 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.

Fig. 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 () followed by the - SN mean RV correlation (). The correlations between the width and the 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.

Figure 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, with the spot being seen face-on at phase 0. The SN mean RV seems to have the smallest spurious variations caused by the spot.
Figure 6: (left) Correlations between the different asymmetry parameters and their corresponding RV estimates in the case of an equatorial 1% spot on the simulated Sun. (right) 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 the 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.

As before, the original RV estimates are corrected using Eq. eq:RV:correction. The results of this correction are displayed in Fig. 7 and the statistical tests on the coefficients involved in Eq. eq:RV:correction 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 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 sdown to 3.04 m s. The corresponding values for SN median RV and SN mean RV are 5.85 m sdown to 2.74  m sand 5.27 m sdown 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. eq:RV:correction are statistically significant to explain the activity signal as seen in Table 3, except the width of the CCF and the interaction term .

Figure 7: (top) The spurious estimated RVs (black dots) caused by a spot in the simulated data, the estimated RVs using Eq. eq:RV:correction (red pluses) and the estimated RVs using the usual correction for stellar activity (green triangles), based on for the SN fit and on for the normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (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 statistical significance on the parameters are presented in Table 3.
Parameter N mean RV SN mean RV SN median RV
Table 3: P–values for the different coefficients used in Eq. eq:RV:correction for the correction of stellar activity induced by an equatorial 1% spot on the simulated Sun. All the parameters corresponding to the Normal or SN parameters are statistically significant to explain the spurious RV variations caused by this spot, except for the width of the CCF and the variable that evaluates the interaction between the width and the asymmetry of the CCF. The proposed correction is not able to perform as well as for the facula, and values for the linear combination are smaller than 0.8.

4.3 Spot and planet

The final simulation includes a planetary signal influencing the CCF along with the 1% spot modelled previously (see Sec. 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 swith no eccentricity and with a period corresponding to one-third of the stellar rotational period, i.e. one-third of 25 days.

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

The correlation between the different CCF parameters are displayed 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 the ones obtained for the spot-only model: –SN median RV has the strongest correlation (), followed by the correlation between BIS SPAN–N mean RV () and finally by the correlation between –SN mean RV (). 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.

Figure 8: RV estimates for N mean RV (red dashed line), SN mean RV (black line) or SN median RV (cyan dotted-dashed line). In this case, the CCFs have been generated using SOAP 2.0, considering an equatorial 1% spot on the simulated Sun in addition to a planet with a period of one-third of the rotational period of the star and with an amplitude of 10  m s. The star does one full rotation between phase -0.5 and 0.5, with the spot being seen face-on at phase 0.
Figure 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.

In order to 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 modelled as a combination of the activity and the planetary signals:

(6)

where can be found in Eq. eq:RV:correction, and , in the case with no eccentricity, can be modelled by the following sinusoidal function:

(7)

with amplitude , orbital period

, and an epoch at the periapsis

. The previous three unknown parameters define the planetary orbit.

The proposed model from Eq. eq:RV:correction.overall was fitted to the RV data and the results of the estimated model are summarized in Table 4. Except for the intercept , the width of the CCF , and the interaction term that evaluates the interaction between the width and the 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 ones (rms = 2.80  m s).

Parameter N mean RV SN mean RV SN median RV
K
P
Table 4: P-values for the different coefficients used in Eq. eq:RV:correction for the correction of stellar activity induced by an equatorial 1% spot on the simulated Sun, and a planet with a period of one-third the rotational period of the star and a semi-amplitude of 10  m s. All the parameters corresponding to the Normal or SN variables are statistically significant, except the intercept , the width of the CCF and the interaction term that evaluates the interaction between the width and the asymmetry of the CCF. Once corrected for stellar activity and the presence of the planet, the residuals are smaller in terms of rms when using the SN variables. Note that because we used nonlinear least squares to fit the best model, the residual rms rather than the is displayed as a reference.

5 Real data application

The analysis 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 Sec. 2.1 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 the Normal and the Skew Normal

A total of CCFs that were derived from the spectra of Alpha Centauri B taken in 2010 by the HARPS spectrograph have been analyzed. 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 signal (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, with a slope equal to and a strong Pearson correlation coefficient of . 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  m s.

Figure 10: (left) Correlation between and the BIS SPAN for Alpha Centauri B. The strong correlation suggests these two parameters are similarly measuring the asymmetry. (top right) RVs as function 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) Differences between the RVs estimated with the Normal density and those from the SN density.

The right plot of Fig. 10 displays 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 Sec. 4, in Fig. 11 we compare the correlations between the asymmetry or the width parameters of the CCF and the RV. For this analysis, we also include the asymmetry parameters derived in Boisse et al. (2011), and in Figueira et al. (2013), BIS-, BIS+, Bi Gauss and , 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 , while the next strongest is 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, with for SN FWHM-SN mean RV, compared to for FWHM-N mean RV.

Figure 11: (top three rows) Correlations between the asymmetry parameters and their corresponding estimated RVs for Alpha Centauri B. (bottom row) The correlation between the FWHM and the estimated RVs. The correlations are stronger when using parameters derived from the SN fit than the Normal one. The estimated ’s are all statistically significant.

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 Sec. 4.3), we observe some significant differences. The correlations between the different parametrisations of the CCF asymmetry and barycentre 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 positives. In the case of the SOAP 2.0 simulated data for a facula or a spot, we always find negatives 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 barycentre 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 modelled 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 barycentre is not possible with current tools, and some progress are still to be made.

Results illustrating the performance of the stellar activity correction proposed in Sec. 3 are displayed 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. eq:RV:correction, 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 Sec. 3 performed only slightly better than the usual correction that uses only a linear combination of the width and the 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 modelling stellar activity. By analyzing the values of the coefficient of determination, , 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 due to faculae. Like observed for the simulated facula in Sec. 4.1, the amplitude of the activity signal is slightly stronger for SN mean RV than for N mean RV or SN 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 Sec. 4.1. Those are arguments strengthen the findings of Dumusque (2014) that also found evidence for faculae dominating the RV stellar signal of Alpha Centauri B.

Figure 12: (top) The RVs (black dots) for Alpha Centauri B estimated using a SN and a Normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (New corr. std, black dots) and the residuals from the usual correction (Usual corr. std, blue triangles), based on for the SN fit and on for the normal fit. The residuals have a smaller systematic component when using the proposed model of Eq. eq:RV:correction (black dots) compared to the usual model (blue triangles).
Parameter N mean RV SN mean RV SN median RV
Table 5: P-values for the different coefficients used in Eq. eq:RV:correction for the correction from stellar activity in Alpha Centauri B data. All the variables corresponding to the Normal or SN parameters are statistically significant, except for the asymmetry of the CCF when using the BIS SPAN with N mean RV. The estimated show that the proposed linear combination explains the most variability in RVs due to the stellar activity when the RVs are estimated with the SN mean RV.

5.2 Comparison for HD192310, HD10700, HD215152 and Corot-7 of the different CCF parameters derived with the Normal and the 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 carryout similar analyses for four other main sequence stars: HD192310 (K2V, Pepe et al. 2011), HD10700 (G8V, Feng et al. 2017), HD215152 (K3V, Delisle et al. 2018) and finally Corot-7 (K0V, Haywood et al. 2014). The same correlation and residual plots displayed 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 HD10700 where the Pearson correlation coefficient between FWHM and N mean RV is equal to , while it is between SN FWHM and SN mean RV, and between SN FWHM and SN median RV. Like in the case of Alpha Centauri B, positive correlations are always observed between the asymmetry and barycentre 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 HD10700, 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, HD10700, HD215152, HD192310 and Corot-7, and the width-RV correlation for Alpha Centauri B, HD215152, HD192310 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 Sec. 3, both RVs estimates had similar residuals. For HD10700, HD215152, and HD192310, 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 the one obtained with the usual correction.

5.3 Detection limits when using the estimated RVs from the Normal or the 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. eq:RV:correction, as shown in the bottom plots of Fig 12, the rms of the residuals are 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, and also compare them when using the usual stellar activity correction with the proposed stellar activity model of Eq. eq:RV:correction. To carryout 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 red-shifted 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 Sec. 4.3, i.e. fitting the activity signal using Eq. eq:RV:correction or the usual correction along with a circular planetary signal (see Eq. eq:RV:correction.overall). Because of the non-linearity of the model that includes a planet, a non-linear 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 here 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 initialised the minimisation 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 was initialised to the time when the RV was crossing 0 since we use a sinusoidal function to fit the planetary signal (see Eq. eq:RV:correction.planet).

Once the parameters involved in Eq. eq:RV:correction.overall were estimated, signals in the residuals, defined as , were analyzed using a Generalized Lomb-Scargle periodogram (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009). If a signal with a P-value555The P-values were estimated using a bootstrap procedure. 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. eq:RV:correction), and when using the usual activity correction. We can see that the new correction for stellar activity based on Eq. eq:RV:correction 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 to detect planetary signals with a slightly smaller amplitude that the usual correction that uses only a linear correlation with the FWHM and BIS SPAN.

Figure 13: Detection limits of planetary signals once the stellar activity signal is removed from the raw RVs using the model proposed in Eq. eq:RV:correction (solid lines) or the usual correction based on for the SN fit and on for the normal fit (dashed lines). The correction for stellar activity based on Eq. eq:RV:correction improves on average the detection limit by 12% and the different RV estimators have similar detection limits.

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 a 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 go to the individual spectrum where each individual points 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 in order 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 HD215152, HD192310 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 the saturation limit of the HARPS detector, respectively. HD10700 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 displayed in Fig. 14. The top row shows the standard errors of the three different estimated RVs, the width, and the asymmetry estimates. Note that because BIS SPAN and do not have the same units, the estimated slopes of the correlation between those two parameters to transform in  m swere used (see Fig. 10 and Table 6 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 will 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 imply 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 use 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.

Figure 14: Results of the bootstrap analyses on the stars HD215152, HD192310 and Corot-7. (top) 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) 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 in average larger than the standard errors of N mean RV (red triangles). However, the standard errors for SN median RV (cyan pluses) are on average smaller than the standard errors coming from the N mean RV. The use of the asymmetry SN parameter leads to standard errors in average smaller than the standard errors related to the BIS SPAN. Note that for asymmetry, the error in BIS SPAN is in  m s. To be able to compare the errors in and BIS SPAN, we multiplied the error in by the slope of the correlation between and BIS SPAN.

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 the 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 the median of the SN density. The width is derived using the variance of the SN density (SN FWHM) 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 the 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 the different SN or Normal parameters. In the case of a spot, the strongest correlations are found between and SN median RV and between BIS 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 between FWHM and N mean RV of R=0.47. The SN parameters continued to have stronger correlations than the Normal ones in the setting where a planetary signal was added the SOAP 2.0 with a single spot.

Looking at real data, we arrive to a similar conclusion that the SN parameters are more sensitive to activity. While real data confirms that the correlation between CCF width and barycentre is always stronger for the SN parameters, they also show that it is the case for the correlations between the CCF asymmetry and barycentre. 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 barycentre. This discrepancy between simulations and real data could be due to the fact that 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 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-barycentre correlation for HD10700, 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 like  Centauri B or HD10700, increasing the photon noise by 60% will not have a significant impact on the RV precision as the instrumental noise will dominate 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, respectively, for the facula and the spot.

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, for the simulations based on the presence of a facula or a spot the new proposed correction almost entirely explains the spurious variations in RV. 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 datasets to see if certain components of the model proposed here 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. eq:RV:correction proposed in this paper to mitigate the effect from stellar activity slightly improves the detection limit with respect to the usual one based on for the SN fit and on 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 precise than N mean RV. Regarding the asymmetry estimates, we observe that has a precision 15% better than the BIS SPAN.

8 Conclusion

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, the width and the 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 that 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 the BIS SPAN. The SN parameters , SN FWHM, and SN mean RV consistently had stronger correlation than those between any of the parameters derived by the Normal and the 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, respectively, 10% and 15% smaller than the corresponding uncertainties on N mean RV and BIS SPAN, though SN mean 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, respectively, for facula and spot. When applying this model on real data, we observe that planetary detection limits are improved by a non-negligible 12%.

9 Acknowledgements

The authors thank Yale’s Center for Research Computing for their help and resources with some of the computational aspects of this work. XD is grateful to The Branco Weiss Fellowship–Society in Science for its financial support. JCK was partially supported by the National Science Foundation under Grant AST 1616086 and by the National Aeronautics and Space Administration under grant 80NSSC18K0443. US 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 Appendix

In this Appendix, a similar analysis as the one presented in Sec. 5 is discussed for four main-sequence stars: HD192310 (K2V, Pepe et al. 2011), HD10700 (G8V, Feng et al. 2017), HD215152 (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 6 summarizes the results obtained by the SN and Normal density models of the CCF. These results 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 or between 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 Fig. 16, 18, 20, and 22, show the corrected RVs. Once corrected from stellar activity, the Normal and SN residuals are comparable for the stars , and . However, the rms of the residuals for Corot-7 are 0.23  m s lower for the SN model than the Normal model. The average S/N at 550 nm for the stars HD10700, HD192310, 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.

Star # CCFs
Star
Table 6: Notable correlations between the asymmetry or the FWHM parameters and the estimated RVs for four stars: , , and . The complete results of the analyses of the correlations for the four stars are presented in Fig. 15, 17, 19, and 21.
Figure 15: (top three rows) Correlations between the asymmetry parameters and their corresponding RVs for . (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 are all statistically significant.
Figure 16: (top) The RVs (black dots) for estimated using a SN and a Normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on for the SN fit and on for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.
Figure 17: (top three rows) Correlations between the asymmetry parameters and their corresponding RVs for . (bottom row) Correlations between the FWHM and the RVs for . 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 the analogous correlations with the estimated SN RVs. The estimated are statistically significant, except for the correlation between FIG BIS and RV (p–values=0.36).
Figure 18: (top) The RVs (black dots) for estimated using a SN and a Normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on for the SN fit and on for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.
Figure 19: (top three rows) Correlations between the asymmetry parameters and their corresponding RVs for . (bottom row) Correlations between the FWHM and the RVs for . The correlations are consistently stronger when using SN mean RV compared to N mean RV. The p–values associated with each 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).
Figure 20: (top) The RVs (black dots) for estimated using a SN and a Normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on for the SN fit and on for the Normal fit. The residuals for both the proposed correction from stellar activity are comparable.
Figure 21: (top three rows) Correlations between the asymmetry parameters and their corresponding RVs for . (bottom row) Correlations between the FWHM and the RVs for . The correlations are consistently stronger when using parameters derived from the SN than the Normal. The p–values associated with each 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).
Figure 22: (top) The RVs (black dots) for estimated using a SN and a Normal fit. (bottom) The residuals from the model fit using Eq. eq:RV:correction (New corr. std–black dots) and the residuals from the usual correction (Usual corr. std–blue triangles), based on for the SN fit and on for the normal fit. The residuals have a smaller systematic component when using the proposed model of Eq. eq:RV:correction (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  m ssmaller than the remaining standard deviation of the Normal model.

References

  • Anglada-Escudé & Butler (2012) Anglada-Escudé, G. & Butler, R. P. 2012, The Astrophysical Journal Supplement Series, 200, 15
  • Arellano-Valle & Azzalini (2008)

    Arellano-Valle, R. B. & Azzalini, A. 2008, Journal of Multivariate Analysis, 99, 1362

  • Azzalini (1985) Azzalini, A. 1985, Scandinavian journal of statistics, 12, 171
  • Azzalini & Capitanio (2014) Azzalini, A. & Capitanio, A. 2014, The skew-normal and related families. Institute of Mathematical Statistics Monographs
  • Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, Astronomy and Astrophysics Supplement Series, 119, 373
  • Belsley (1991) Belsley, D. A. 1991, Conditioning diagnostics: Collinearity and weak data in regression No. 519.536 B452 (Wiley New York)
  • Boisse et al. (2011) Boisse, I., Bouchy, F., Hébrard, G., et al. 2011, Astronomy & Astrophysics, 528, A4
  • Boisse et al. (2009) Boisse, I., Moutou, C., Vidal-Madjar, A., et al. 2009, Astronomy & Astrophysics, 495, 959
  • Borgniet et al. (2015) Borgniet, S., Meunier, N., & Lagrange, A.-M. 2015, A&A, 581, A133
  • Bouchy et al. (2005) Bouchy, F., Pont, F., Melo, C., et al. 2005, A&A, 431, 1105
  • Cavallini et al. (1985) Cavallini, F., Ceppatelli, G., & Righini, A. 1985, Astronomy and Astrophysics, 143, 116
  • Claret & Bloemen (2011) Claret, A. & Bloemen, S. 2011, A&A, 529, A75
  • Cosentino et al. (2012) Cosentino, R., Lovis, C., Pepe, F., et al. 2012, in Proc. SPIE, Vol. 8446, 84461V
  • Delisle et al. (2018) Delisle, J.-B., Ségransan, D., Dumusque, X., et al. 2018, A&A, 614, A133
  • Desort et al. (2007) Desort, M., Lagrange, A.-M., Galland, F., Udry, S., & Mayor, M. 2007, Astronomy & Astrophysics, 473, 983
  • Dumusque (2014) Dumusque, X. 2014, ApJ, 796, 133
  • Dumusque (2016) Dumusque, X. 2016, Astronomy & Astrophysics, 593, A5
  • Dumusque (2018) Dumusque, X. 2018, ArXiv e-prints [[arXiv]1809.01548]
  • Dumusque et al. (2014) Dumusque, X., Boisse, I., & Santos, N. 2014, The Astrophysical Journal, 796, 132
  • Dumusque et al. (2017) Dumusque, X., Borsa, F., Damasso, M., et al. 2017, Astronomy & Astrophysics, 598, A133
  • Dumusque et al. (2012) Dumusque, X., Pepe, F., Lovis, C., et al. 2012, Nature, 491, 207
  • Dumusque et al. (2011) Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. 2011, Astronomy & Astrophysics, 525, A140
  • Faria et al. (2016) Faria, J., Haywood, R., Brewer, B., et al. 2016, Astronomy & Astrophysics, 588, A31
  • Feng et al. (2017) Feng, F., Tuomi, M., & Jones, H. R. 2017, Astronomy & Astrophysics, 605, A103
  • Feng et al. (2017) Feng, F., Tuomi, M., Jones, H. R. A., et al. 2017, AJ, 154, 135
  • Figueira et al. (2013) Figueira, P., Santos, N., Pepe, F., Lovis, C., & Nardetto, N. 2013, Astronomy & Astrophysics, 557, A93
  • Fischer et al. (2016) Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, Publications of the Astronomical Society of the Pacific, 128, 066001
  • Hatzes (2002) Hatzes, A. P. 2002, Astronomische Nachrichten, 323, 392
  • Haywood et al. (2014) Haywood, R., Collier Cameron, A., Queloz, D., et al. 2014, Monthly notices of the royal astronomical society, 443, 2517
  • Herrero et al. (2016) Herrero, E., Ribas, I., Jordi, C., et al. 2016, A&A, 586, A131
  • Kurster et al. (2003) Kurster, M., Endl, M., Rouesnel, F., et al. 2003, ASTRONOMY AND ASTROPHYSICS-BERLIN-, 403, 1077
  • Lagrange et al. (2010) Lagrange, A.-M., Desort, M., & Meunier, N. 2010, Astronomy & Astrophysics, 512, A38
  • Levenberg (1944) Levenberg, K. 1944, Quarterly of applied mathematics, 2, 164
  • Lindegren & Dravins (2003) Lindegren, L. & Dravins, D. 2003, Astronomy & Astrophysics, 401, 1185
  • Lomb (1976) Lomb, N. R. 1976, Astrophysics and Space Science, 39, 447
  • Marquardt (1963) Marquardt, D. W. 1963, Journal of the society for Industrial and Applied Mathematics, 11, 431
  • Mayor et al. (2003) Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20
  • Meunier et al. (2010) Meunier, N., Desort, M., & Lagrange, A.-M. 2010, Astronomy & Astrophysics, 512, A39
  • Oshagh et al. (2013) Oshagh, M., Boisse, I., Boué, G., et al. 2013, A&A, 549, A35
  • Pepe et al. (2011) Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58
  • Pepe et al. (2002) Pepe, F., Mayor, M., Galland, F., et al. 2002, Astronomy & Astrophysics, 388, 632
  • Pepe et al. (2014) Pepe, F., Molaro, P., Cristiani, S., et al. 2014, Astronomische Nachrichten, 335, 8
  • Queloz et al. (2009) Queloz, D., Bouchy, F., Moutou, C., et al. 2009, Astronomy & Astrophysics, 506, 303
  • Queloz et al. (2001) Queloz, D., Henry, G., Sivan, J., et al. 2001, Astronomy & Astrophysics, 379, 279
  • Rajpaul et al. (2015) Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, Monthly Notices of the Royal Astronomical Society, 452, 2269
  • Robertson et al. (2014) Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, Science, 1253253
  • Saar & Donahue (1997) Saar, S. H. & Donahue, R. A. 1997, The Astrophysical Journal, 485, 319
  • Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  • Teunissen (1990) Teunissen, P. 1990
  • Thompson et al. (2017) Thompson, A., Watson, C., de Mooij, E., & Jess, D. 2017, Monthly Notices of the Royal Astronomical Society: Letters, 468, L16
  • Zechmeister & Kürster (2009) Zechmeister, M. & Kürster, M. 2009, A&A, 496, 577