© N. Meunier et al. 2019

1 Introduction

The detection of low-mass planets using the radial velocity (RV) technique is strongly impacted by the presence of stellar variability. Magnetic activity leads to spurious RV signals around the rotational period as well as on longer timescales (related to cycle variations). We have shown that on long timescales the signal is mostly due to the inhibition of the convective blueshift in plages (Meunier et al. 2010a), producing not only a strong peak in the periodograms at the cycle period, but also a significant power at higher frequencies (periods of a few hundred days and above): this decreases the detection limit of low-mass planets by one to two orders of magnitudes. A commonly used technique to correct for this contribution relies on the strong correlation between this RV component and the chromospheric emission as measured by the classical , since both depend on the filling factor covered by plages. Such an approach has been used in a large number of publications (e.g. Boisse et al. 2009; Pont et al. 2011; Dumusque et al. 2012; Robertson & Mahadevan 2014; Rajpaul et al. 2015; Lanza et al. 2016; Díaz et al. 2016; Borgniet et al. 2017) and allows a significant decrease in the RV jitter to be achieved. This correlation is also used in the fitting challenge organised by X. Dumusque to compare the performance of up-to-date correction techniques (Dumusque 2016; Dumusque et al. 2017).

The limitations of this approach are however not well studied. In the solar case, we have shown that although this approach allows a significant gain, it is not suitable for very low planet masses (Meunier & Lagrange 2013): for example, at distances of around 1 AU, planets of 1 M_Earth could only be reached with very good temporal samplings. It is critical to understand how these limitations depend on the star (spectral type, activity level, ...) and why a better performance cannot be achieved with this method. The issues are the following: (1) the dependence on can be fitted to a certain extent by a polynomial in time, which is often used to remove the contribution from possible unknown companions at very long periods that lead to a degeneracy between the and this polynomial. (2) The correction using the time series is not perfect, and a significant amount of power remains in the periodograms, especially at long periods (around the cycle period and below).

In this paper, we investigate these issues using a large set of simulations covering F6–K4 old main sequence stars with various activity levels. The generation of the time series is described in Meunier et al. (2019), hereafter referred to as Paper I, which also shows how the , representingthe chromospheric emission (average level, cycle amplitude) varies with inclination, with results which agree with Shapiro et al. (2014). An initial analysis of the RV jitter was made in Meunier & Lagrange (2019a), who compared the RV jitter from simulations with observed jitter (Saar & Brandenburg 1999; Santos et al. 2000; Wright 2005; Isaacson & Fischer 2010). They also showed the good agreement with the RV versus spectral type slope observed by Lovis et al. (2011) on a large sample of HARPS observations. Finally, they showed that the current status of correction techniques would not allow the detection of Earth-mass planets around stars like the Sun, and that for lower-mass stars a very large number of points was necessary to reach that goal. A significant improvement of the correction method must therefore be made. Physical models are needed to remove the activity contribution to ensure that the residuals are well controlled: the correction to be made is typically one to two orders of magnitude (depending on the activity level). We show in the present paper that the time series based on complex and realistic activity patterns can help us to find clues to develop better models to perform these corrections, which in the end should help to better control the residuals after correction.

The outline of the paper is the following. In Sect. 2 we briefly recall the model and parameters. In Sect. 3, we study the correlation between RV and time series across our grid of parameters. Subsequently, we present a study of the typical gain which is obtained when performing standard correction using and a polynomial in time. In Sect. 5, we identify a limiting element in the correction and characterize it. We propose a new method to correct for this effect in Sect. 6, which will impact the detectability at long orbital periods. Finally, we conclude in Sect. 7.

2 Model and parameters

Our model, described in detail in Borgniet et al. (2015) for the Sun and in Paper I for other stars (F6–K4 and relatively old main sequence stars), provides consistent spots, plages, and network structures for complex activity patterns similar to the Sun.

The time series due to activity is the sum of the contributions due to spots (rvspot1 and rvspot2, where two laws are used for the spot temperature contrast¹), plages rvplage (contrast dependence on spectral type from Norris 2018), and inhibition of the convective blueshift in plages rvconv. We also consider the addition of the contribution of oscillation, granulation, and supergranulation (OGS) signal, rvogs, averaged over 1 h, and an instrumental white noise represented by a Gaussian noise with an amplitude of 0.6 m s⁻¹ (rvogs). Here, time series are produced to be able to relate RV variation to chromospheric emission. We also produced photometric (Meunier & Lagrange 2019b) and astrometric (Meunier & Lagrange, in prep.) time series, but the analysis of these time series is outside the scope of the present paper as we focus on the relationship between RV and chromospheric emission. The temporal step is one day on average (with random departures of up to 4 h). The time series have a maximum length of 15 yr, and always cover an integer number of cycles.

The rotation periods, and the cycle periods and amplitudes also depend on the star, and a range of realistic values is considered for each spectral type and activity level. Some parameters which are not constrained are kept to the solar values considered in Borgniet et al. (2015), for example meridional circulation or size distributions. We refer to Paper I for more details about the laws used to produce the consistent sets of parameters, where all references and justifications can be found. Of particular interest in this paper is the maximum average latitude at the beginning of cycle θ_max: it is not yet constrained from observations or models, and three values were considered in Paper I, the solar latitude θ_max,⊙, θ_max,⊙+10°, and θ_max,⊙+20°, so that the activity pattern covers different ranges in latitude.

3 Correlation between RV and log R

In the solar case, we showed that inhibition of the convective blueshift was dominating over the spot and plage contributions (Meunier et al. 2010a; Borgniet et al. 2015). We know from other stars that the correlation is not always strong, as seen for example from the RV slope in Lovis et al. (2011): the slope can be small in some cases (see also Meunier et al. 2017), and it is important to better understand and quantify how this property varies with spectral type or inclination for example to interpret the observations. It is usually assumed that departures from a perfect correlation are due to the addition of the spot and plage signal for example, but this may not be the only cause.

We first consider the correlation² between RV and , because it can easily be derived from observational data and gives some clues on what to expect from the correction. Typical values are shown in detail in Appendix A.1. The correlation between RV and depends on many parameters, but the most influential are spectral type, activity level, and inclination. The departure from a correlation of one explains why using to correct RV times series is not perfect. In addition, the presence of noise (OGS or instrumental noise) strongly degrades the global correlation. It is for example possible to have a RV jitter dominated by rvconv, although the correlation reaches relatively low values (Appendix A.2). This is probably due to the fact that in some cases, the addition of the OGS signal significantly degrades the correlation even though its contribution to the RV jitter is not major.

The global correlation between RV and includes contributions from both short and long timescales. The short-term correlations presented in Appendix A.3are related to the global correlation, but there are a large number of simulations with short-term correlations that are much lower than the global one; these are always lower than 0.75. The inhibition of the convective blueshift is less dominant even for edge-on stars, and there are many configurations where the inhibition of the convective blueshift is not dominant on short timescales. In addition, weobserve a very strong effect of inclination: for pole-on stars, the short-term correlations are much lower than the global ones. This is due to the fact that for low inclinations the rotation modulation is much smaller and therefore the short timescale signal becomes noisier.

We also note that if the inclination of a given star is well known, it seems possible to relate the global and average short-term correlation. The observation of a given short-term correlation at a given time is however not representative of the average short-term correlations, and its sign can even be negative as local correlations cover a wide range of values. Finally, local correlations computed on solar time series also show a large dispersion similar to what is observed in our simulations (Appendix A.4).

	Fig. 1 Gain in RV jitter (when correcting for alone) vs. global correlation between RV and (including rvogs and rvinst).
Open with DEXTER

4 Decrease of the gain in RV jitter when correcting using the log R–RV correlation

To characterize the gain in RV jitter, we use the ratio between the RV jitter before and after correction. The ratio should be higher than one to get an efficient correction, and the higher the ratio, the better the correction. We first consider the relationship between gain and global correlation, and analyse how the gain varies for different types of corrections. The first model used to perform the correction is of the form RV. In a second step, we discuss the addition of a polynomial in time, described in the introduction: the model is then of the form RV (Dumusque et al. 2017).

4.1 Gain while correcting from log R and polynomial in time

It is interesting to see whether it is easier to correct for activity using the –RV correlation when the correlation between RV and is high. Figure 1 shows the gain defined at the beginning of Sect. 4 versus the –RV correlation (from Sect. 3) and indeed shows a perfect relationship between the two approaches. When the correlation departs from one, the gain also decreases.

The gain as a function of the original RV jitter (before correction) is shown in Fig. 2. The upper panel corresponds to a simple correction using a linear relation between RV and and we see a complex behaviour with a wide spread: a large number of simulations correspond to low to medium RV jitter and a gain close to the solid line (gain providing a RV jitter of below 1 m s⁻¹ after correction). For RV jitters in the range 1–2.5 m s⁻¹, 37% can be corrected to a level below 1 m s⁻¹. However, for medium to high RV jitters, the gain can be very different for a given RV jitter: for example, for a jitter of 4 m s⁻¹, the gain can vary between 1 (no improvement brought by the correction) and 3.5. The few simulations with the highest RV jitter are very poorly corrected (gain close to 1). Simulations with a RV jitter originally above 2.5 m s⁻¹ never reach RV jitter below 1 m s⁻¹ after correction. Finally, the shape of the upper envelope shows that the gain does not vary linearly with RV jitter.

The second panel shows the same plot, but this time the RV signal is corrected for both and a second degree polynomial in time (which is often doneto correct for long-term trends due to companions). The general behaviour is similar, although it is possible in some cases to obtain higher gains, as illustrated on the last panel, which shows the gain with and polynomial in time correction vs. gain with correction only. This is puzzling, because in the present case, the polynomial in time is not a good model of the time series, since it is not present in our simulations (no added companions): when correcting using this polynomial in time in addition to the linear correlation with , we use a model which is not physical.

Fig. 2 Gain in RV jitter in various conditions. Upper panel: gain vs. RV jitter (for simulations including rvogs and rvinst, and for ΔTspot2), for a correction with alone. The straight black line has a slope of 1, indicating points where the correction allows a residual RV jitter of 1 m s−1 to be reached. The colour code corresponds to inclination, from pole-on (i = 0°, yellow) to edge-on (i = 90°, blue), with light and dark orange corresponding to 20° and 30°, light and dark red to 40° and 50°, brown to 60°, and light and dark green to 70° and 80°. Middle panel: same for a correction with and a second-degree polynomial in time. Lower panel: gain with and time correction vs. gain with correction only. Only one point out of five is shown for clarity.

Open with DEXTER

	Fig. 3 Upper panels: amplitude of the polynomial in time vs. B−V (left) and (right). The black dots represent each individual time series, while the red points correspond to a binning in B− V and respectively. Lower panels: percentage of polynomial amplitudes higher than 0.5 m s−1 (stars), 1 m s−1 (diamonds), and 2 m s−1 (triangles) vs. B−V and respectively.
Open with DEXTER

4.2 Discussion on the polynomial function fitted to the RV time series

Given the difference in gain observed for certain simulations when correcting for alone and when also introducing a polynomial in time, it is important to quantify this effect and understand its origin, especially since it is widely used. Furthermore, the simulations do not include companions and therefore the coefficients of the polynomial should be equal to zero (i.e. the polynomial vs. time should be flat); if the coefficients of the polynomial in time are different from zero, this means that we fit the time series with a model which is not physically realistic.

For each simulation, we therefore corrected for and this polynomial in time as above. We then computed the amplitude of the fitted polynomial A_pol, defined as the maximum minus the minimum of the polynomial over the time range. Figure 3 shows this amplitude as a function of B− V and (upper panels). Although the average is low, there are simulations with high values of A_pol, up to a few m s⁻¹. The percentage of simulations with amplitudes above 0.5, 1, and 2 m s⁻¹ is shown in the lower panels. For the higher masses and more active (from their average activity level) stars in the lower panels of Fig. 3, the percentage can be relatively high, for example more than 30% for the 0.5 m s⁻¹ threshold. The percentage can reach a few percent for the 2 m s⁻¹ threshold.

All simulations with high polynomial amplitudes correspond to long cycles, that is, simulations for which only one cycle is simulated, and the highest amplitudes correspond to high cycle amplitudes. For such simulations, the polynomial in time can mimic activity (at least to a certain point, since the shape of the cycles is not exactly polynomial). Since some fits are better with this polynomial than without, the RV variations contain a component which is closer to such a polynomial (as far as a single cycle is concerned) than the variability: a polynomial fit is then ad hoc and the time series are not well constrained in these cases. We explain this result in the following section.

5 Hysteresis between RV and log R time series

In this section, we explain the previous results by the presence of an hysteresis pattern between the RV and chromospheric emission variabilities. We quantify the amplitude of this pattern as a function of the stellar parameters in our simulations. Finally, we show that it is also observed for the Sun and other stars.

5.1 Why does the polynomial in time improve the correction?

We first show an example to illustrate the potential impact of the polynomial fit. All plots are for ΔT_spot2, and includeactivity, OGS (smoothed over 1 h), and instrumental white noise. Figure 4 shows an example of a time series for a G2 star with a medium activity level: all points (full cycle) and after smoothing for clarity. The upper panel shows the difference between the RV we want to correct and the RV that would be obtained if a linear relationship between RV and were used. The lower panel shows RV versus : although there is indeed a very good global correlation between the two (around 0.9 in this example, hence the usual correction method), there is a departure from a strict correlation, with an asymmetry between the ascending phase of the cycle and the descending phase, particularly visible for the pole-on (right panels) plot. For the star seen edge-on, for a given level, the RV is higher during the descending phase of the cycle compared to the ascending phase. It is the opposite for the same star seen pole-on. This is a major effect because we observe differences higher than 1 m s⁻¹ in these examples, which should lead to an important power in the periodograms of the time series.

There is therefore a kind of hysteresis (we use this denomination in the following), with the two phases of the cycle behaving in a different manner. This effect is not taken into account when modelling RV as a linear function of , but in some cases can be partially taken into account using the polynomial in time. There also might be a systematic effect of inclination. We characterize this hysteresis in the following section for all simulations.

5.2 Characterisation of the hysteresis across the grid of parameters

In order to characterise the hysteresis, we define the criterium C_hyst. For a given simulation, we consider one cycle. We estimate the position of cycle maximum to separate the ascending phase and the descending phase. We then consider a series of 30 equally spaced levels (corresponding to the range covered by each time series). The RV is averaged over each bin during the ascending phase (RV_asc) and during the descending phase (RV_desc). We then compute RV_desc-RV_asc and average it over the 30 levels to produce C_hyst.

Figure 5 shows C_hyst versus B−V, , and RV jitter. A first striking result is that there is indeed a systematic effect of inclination, going from a negative (pole-on) to a positive (edge-on) value. The reversal happens around 60° from pole-on. C_hyst is strongly related to the average activity level and RV jitter (i.e. variability level). Because of this, C_hyst naturally decreases towards lower masses.

The value of θ_max (maximum average latitude at the begining of the cycle) also has a strong effect on the hysteresis. Figure 6 shows binned C_hyst vs. the RV jitter, for our three values of θ_max. In addition to an almost linear relationship with the RV jitter in most cases, θ_max strongly impacts C_hyst as well: higher θ_max means a wider range covered by C_hyst. Finally, the amplitudes are lower for edge-on configurations compared to pole-on configurations.

In summary, C_hyst is strongly related to the activity level, inclinations, and θ_max. We attribute this dependence to the conjunction of two facts: structures are not at the same position in latitude on the disk during the cycle, and projection effects are different for RV and . This is detailed and discussed in Sect. 6, where a new correction method is proposed.

	Fig. 4 Example of smoothed time series for a moderately active G2 star seen edge-on (left) and pole-on (right). First line: RV vs. time (solid line), and RV fitted from the variations (dashed line). Second line: vs. time. Third line: RV vs. , during the ascending phase of the cycle (red) and during the descending phase (green).
Open with DEXTER

5.3 Is the hysteresis present in observations?

We show that the hysteresis is also observed for the Sun, based on the analysis of two cycle-long time series. We then find that it is also present for other stars.

5.3.1 The solar case

Our simulations provide evidence for the presence of an hysteresis pattern between the long-term RV and variations over the cycle. We now examine two solar RV time series to confirm this property with observations.

We first consider the RV reconstruction over a complete cycle made by Meunier et al. (2010a). This reconstruction was based on observed solar structures (spots, plages), and a model was used to build the integrated RV from the estimated RV for each structure. After selecting days for which an observation of the chromospheric emission was available (S-index from the Sacramento Peak Observatory), we plot the hysteresis pattern in Fig. 7 (upper panels). The amplitude of the hysteresis is of the order of 0.5 m s⁻¹, which is very similar to our simulation for edge-on configurations, both in sign and amplitude.

We also computed the hysteresis pattern from the solar RV time series reconstructed from SOHO/MDI Dopplergrams (Scherrer et al. 1995) by Meunier et al. (2010b). In this case, the RV due to active regions was reconstructed by integrating the Doppler velocities over the disk; we did not compute the RV using the model for the RV associated to each structure. The chromospheric emission is also from the Sacramento Peak Observatory. The hysteresis shown in Fig. 7 (lower panels) is very similar to the previous case, which shows that the model is very good to that level of detail.

These two solar series are, to our knowledge, the only ones available to compare with our simulation. The ongoing solar programs with HARPS-N (Dumusque et al. 2015; Collier Cameron et al. 2019; Milbourne et al. 2019) and HARPS produce high-cadence observations of the Sun in stellar conditions, and will be suitable for such an analysis, but the temporal coverage so far is still insufficient to allow such a study (only the end of the descending phase of the solar cycle is available, with a solar minimum in2019). The reconstruction of the solar RV from structures derived from HMI observations by Milbourne et al. (2019) is similar in length.

5.3.2 Stellar observations

It is more difficult to check whether such an hysteresis is observed on other stars because of the usually poor temporal sampling. Surveys such as the Mount Wilson survey provide large samples of stars with a good cycle coverage (e.g. Baliunas et al. 1995), but they are not associated to simultaneous RV measurements. We have however identified eleven stars observed with HARPS for which a complete cycle with a large number of points is available (Table 1). Unlike in our simulations, the cycle is not necessarily covered from one minimum to the next. We binned the data over one year, and plot the corresponding hysteresis pattern. The results are shown in Fig. 8. We observe both types (edge-on or pole-on sign) of hysteresis in seven stars, and four stars show either a more complex mixed pattern (two stars) or no hysteresis (two stars). This is summarised in Table 1. The presence of stars with no hysteresis is expected given the reversal observed in Fig. 5. A mixed pattern could be due to a spatio-temporal distribution of the structures that is more complex than the solar one. We therefore find evidence suggesting that the same hysteresis is present in other stars as well.

	Fig. 5 Upper panels: hysteresis amplitude Chyst vs. B− V (left), vs. (middle), and vs. RV jitter (right), for F stars. The colour code is similar to that in Fig. 2. Middle panels: same for G stars. Lower panels: same for K stars. Only one point out of five is shown for clarity.
Open with DEXTER

6 Towards a better correction of the long-term RV variability

In this section, we explain the origin of the hysteresis pattern. We then use these results to propose a new correction method using a better model for the relationship between RV and , and illustrate itsperformance on a subset of G2 star simulations.

6.1 Explanation of the inclination-dependent hysteresis pattern

Our interpretation of this hysteresis pattern is that although RV and chromospheric emission are roughly correlated in the long term due to the major contribution of rvconv, the signal produced by a given active region strongly depends on its position on the disk, while the projection effects are different for both variables: the chromospheric emission is essentially linear vs. the projected area of the structure. In RV, the dependence on μ (cosine of the angle between the normal to the surface and the line-of-sight) is more complex, as it includes an additional projection effect of the velocity field (μ), a contrast dependence on μ, and the center-to-limb darkening function (see Paper I).

As a consequence, if the average μ (hereafter ) changes in time (in particular due to a butterfly diagram pattern), this must introduce a departure from the linear relationship between RV and chromospheric emission. This is the case when considering different stellar inclinations. For a solar-like butterfly diagram such as in our simulations (i.e. toward the equator), if the star is seen pole-on, decreases during the cycle because the structures appear at lower latitudes, as shown in Fig. 9, and therefore also covers a wide range: RV decreases faster than the chromospheric emission, so for a given activity level, the RV signal is higher during the ascending phase of the cycle. This is the opposite for a star seen edge-on because during the cycle, increases and covers asmall range of values: for a given activity level, the RV signal is lower during the ascending phase of the cycle. This behaviour should be responsible for the reversal seen in Fig. 5 and is related to the average position of the structures, which is due to the presence of a butterfly diagram pattern, and therefore to dynamo processes. It could also explain why θ_max has a strongeffect on the hysteresis pattern.

Furthermore, a simulation similar to the one shown in Fig. 4, but with a constant latitude over time (flat butterfly diagram), exhibits no hysteresis, as illustrated in Fig. 10. This shows that the hysteresis is strongly related to the spatio-temporal distribution of the structures. We also note from this figure that the relationship between RV and , although showing little dispersion, is not strictly linear and some dispersion is still present.

	Fig. 6 Average hysteresis vs. RV jitter for three values of θmax: solar (black), solar +10° (red), and solar +20° (green). Each plot corresponds to a different inclination, from pole-on (0°, upper left panel) to edge-on (90°, lower right panel). All plots are on the same scale in m s−1.
Open with DEXTER

Fig. 7 Relative velocity (solid line), RV deduced from the linear relation with Sacramento Peak S-index (dashed line) vs. time (left panels), and RV vs. S-index (right panels, in red for ascending phase and in green for descending phase). Upper panels: RVs reconstructed from observed structures (Meunier et al. 2010a). Lower panels: RVs reconstructed from MDI/SOHODopplergrams (Meunier et al. 2010b).

Open with DEXTER

	Fig. 8 Relative velocity vs. time (stars linked by solid black lines), and RV from linear dependence on (diamonds linked by orange dashed line) for each star. The second plot for each star shows the hysteresis pattern, with the ascending phase in red and descending phase in green.
Open with DEXTER

	Fig. 9 Average vs. cycle phase for pole-on (solid line) and edge-on (dashed line) configurations for simulation shown in Fig. 4.
Open with DEXTER

	Fig. 10 Relative velocity vs. for simulation similar to Fig. 4, seen pole-on (upper panel) and edge-on (lower panel), but with constant latitude over time.
Open with DEXTER

6.2 Construction of a reference catalogue

If the butterfly diagram of a star were known, the difference in trend between RV and chromospheric emission could then be modelled and corrected. In practice however, it is not known, but we can build a large number of functions corresponding to different configurations (different latitude ranges for example), which can be generated to build a reference catalogue describing the different possibilities for the ratio between the RV and behaviour over time. This must then be modulated by the activity level over time when applied to a given time series. The function leading to the lowest residuals after correction is then selected. We built these functions for the following parameters:

• Stellar inclination: we consider 91 values between 0° and 90° with a step of 1°.

• Average latitude at the beginning of the cycle θ_max: we consider values between 15° and 59°. We recall that in our simulations, we only considered input values of 22°, 32°, and 42° (the actual average latitudes are usually higher, in particular due to meridional circulation). Here we generate the functions for a wider range of parameters since θ_max is not constrained for a particular observation.

• Width of the butterfly diagram in latitude, Γ: we considered values from 2 to 20°, with a step of 2°. We recall that all simulations were made with an input value of 6° (the actual values are usually higher due to diffusion and meridional circulation). We use a Gaussian distribution around the average latitude, with a cut at ±20°.

• Average latitude at the end of the cycle θ_min: this is kept constant (9°, corresponding to the input value in our simulations) in this first analysis.

For each of these parameter sets, we computed the butterfly diagram pattern in latitude for 100 phases during a cycle. We attributed to each pixel of the stellar disk the projection effects corresponding to RV and chromospheric emission respectively. The center-to-limb darkening is similar to the function used in our simulation in Paper I (Claret & Hauschildt 2003). The plage contrast is an average as a function of μ of the function used in Paper I. For a given inclination, adding the pixel contributions over the whole disk at each phase produces two time series, hereafter RV_cat and Ca_cat. These functions do not include any structures such as spots or plages: they only describe the relative variability of the two variables (RV and chromospheric emission), for a constant activity level, due to the position of the structure. For that reason, we also attribute to each pixel a factor describing the fact that a structure rotating in longitude would spend more time at a position close to the limb compared to disk center due to projection effects (for a given rotation rate). This leads to a catalogue of 20 930 functions. Figure 11 shows an example of such a function (upper panels).

6.3 Fitting the time series using the reference catalogue

The catalogue is then used as follows for a given time series covering a stellar cycle:

• We first bin the time series over the rotation period, because the catalogue is used to improve the correction over long timescales (we average the structure positions in longitude). This leads to RV_bin and Ca_bin vs. t_bin in the following. In principle, the same procedure could be applied directly to the original time series, but this approach is faster and gives a very good idea of the performance when considering the long-term variations only.

• Each function of the catalogue is then tested: we apply a linear fit to (RV_bin, Ca_bin), leading to a RV time series corrected from the linear correlation with Ca_bin, RV_corlin. Here, RV_corlin is then multiplied by RV_cat/Ca_cat (after interpolating on t_bin) to account for the difference in behaviour with μ over the cycle of RV and chromospheric emission, and by a constant to minimize the residuals, the amplitude in the catalogue being arbitrary.

• Finally, we select the catalogue function providing the lowest rms of the residuals.

When applied to the binned time series corresponding to Fig. 4, the rms RV before correction is 3.79 m s⁻¹ (3.65 m s⁻¹), after a standard linear correction 0.48 m s⁻¹ (0.77 m s⁻¹), and after the new correction 0.2 m s⁻¹ (0.27 m s⁻¹), for the edge-on and the pole-on configurations respectively. There is therefore a gain of 2.40 (2.85) between the standard correction and the new correction.

To illustrate the performance of the new method, we applied this procedure for a subset of G2 star simulations, over one cycle only and ΔT_spot1. We compare the rms before correction (rms₀), after a standard linear correction with chromospheric emission (rms_stand), and after this new correction method (rms_new). Median values of the gain are shown in Table 2 and distributions in Fig. 12. With the new method, the gain with respect to the standard correction is about 2.5 for pole-on configurations and 1.5 for edge-on configurations (median values), and very close to one for inclinations of 60° as shown in panel C. This is expected because the reversal of the hysteresis occurs around this inclination, meaning that no improvement is expected in that case. The distributions show that for some simulations the gain is as high as 4–6 with respect to the standard correction.

Fig. 11 First panel: example of correcting function from catalogue (see Sect. 6.4) vs. phase for edge-on configuration (left) and pole-on configuration (right): chromospheric emission (solid line), RV (dashed line), and ratio between RV and chromospheric (thick dotted-dashed line), all normalized to 1 for clarity. Second panel: RV vs. chromospheric emission for the simulation shown in Fig. 4, colour-coded according to the cycle phase: red for the ascending phase, green for the descending phase. The orange and brown curves correspond to the same phases respectively but for the model derived from the catalogue functions. The black dashed line is the linear fit between RV and . Third panel: RV residual after standard correction (black dashed line) and after the new correction (orange and brown solid lines as in the second panel). Fourth panel: RV residual vs. time after standard correction (dashed line) and after the new correction (solid line).

Open with DEXTER

	Fig. 12 Distribution of different rms ratios, for inclinations of 0° (orange), 30° (red), 60° (brown), and 90° (green). Panel A: rms0/rmsstand. Panel B: rms0/rmsnew. Panel C: rmsstand/rmsnew. Panel D: rmsnew/rmsopt. Panel E: rmsnew/rmsext. Panel F: rmsext/rmsopt.
Open with DEXTER

6.4 Going beyond the hysteresis correction

The residuals show that no hysteresis remains, but a curvature is still present (Fig. 11). This is likely to be related to the non-linearity observed in Fig. 10. After correction of this curvature using a second-degree polynomial in on the residuals shown in Fig. 11, the residuals have a rms of 0.219 and 0.174 m s⁻¹ (for pole-on and edge-on respectively): this is very close to the residuals after a similar correction (second degree polynomial) for the time series at constant latitude shown in Fig. 10 (rms of the residuals of 0.215 and 0.176 m s⁻¹). We propose to add another step to the procedure, consisting in performing this second degree in the fit, which leads to an rms of the residuals, rms_ext. The results are shown in Fig. 12 (panels E and F) and Table 2: for many simulations, the gain is centred on one, but it allows a significant gain to be added for some of them (tail in panel E).

After this step, the median gain in power for periods in the habitable zone of G2 stars computed as in Meunier & Lagrange (2019a), that is in the 274-777 days range, is between 1.3 (60°) and 2.8 (30°). The gain is higher than two (representing a gain of four in mass) for 61, 63, 21, and 49% of the simulations for the four inclinations respectively. We note that for 7–18% of the simulations, the power is slightly increased after correction however (the gain is very high for periods between typically P_cyc/2 and P_cyc, but not as high below P_cyc).

As a complementary approach, we also fitted the time series with a function of the form RV =α(1 + + γ()²)(1 + δt + ϵt²) where the polynomial in takes a non-linear relationship between RV and chromospheric emission into account, and the polynomial in time plays the role of the reference catalogue functions computed above to correct for the hysteresis. This represents an estimate of the best gain which can be achieved when considering processes on long (cycle) time scales, as it is less physically constrained than the current functions in the catalogue. We computed the rms of the residuals after such a correction, rms_opt. The distribution of the gains are shown in Fig. 12 (panels D and F) and the median gains in Table 2. The new correction (previous section) provides a gain with respect to rms_opt with a peak around 1, and a small tail toward higher values. After the polynomial fit (rms_ext), the tail has disappeared, which shows that the correction is very close to what we can expect to reach, which is an excellent result.

We note that at this stage, the rms of the binned series after correction is in the range 0.02–0.33 m s⁻¹. There is therefore some signal left, which is more stochastic in nature than what we have removed so far, and is on typically lower timescales. We will study these residuals in a future work. They could be due to departure from the average latitude of large active regions, and/or to size effects (see Fig. 10 and Sect. 6.1).

In conclusion, this section shows the feasibility and interest of this method, which appears to be very promising and allows substantial improvement of the residuals. In a future work, we will test the impact of a degraded sampling (lower number of points, incomplete coverage of the cycle). The construction of the catalogue will also be improved using smaller steps and/or interpolating in the catalogue, and taking different values of θ_min into account. Finally, the performance will also be tested when a planet in the habitable zone is added, and for all spectral types in our simulations (F6–K4).

7 Conclusion

We analysed a large number of simulated time series of stellar activity covering the spectral types F6–K4 and different activity levels. We showed that such simulations, based on realistic complex activity patterns, are very useful to understand the limitation of correction methods. We have focused on the methods using the correlation between RV and .

The detailed study of the correlation between RV and , and the gains in RV jitter which can be obtained using the typical correction compared to the gain obtained when also considering a second-degree polynomial in time, led us to several important conclusions:

• Inclination usually plays a crucial role.

• The global and short-term correlations related to the relative weight between rvspot and rvplage (contribution of spots and plages due to their intensity contrast) on one side, and rvconv (due to the inhibition of the convective blueshift in plages) on the other side, present a complex relationship. This relation is impacted by inclination, and by the addition of other contributions such as oscillation, granulation, and supergranulation.

• Not only is the correction using limited (as shown in previous papers Meunier & Lagrange 2013, 2019a; Dumusque et al. 2017), but the addition of a polynomial in time may also help to reduce the RV jitter a little, although in an uncontrolled way. However, the results show that adding polynomial fits in time should be used with caution.

• An hysteresis between RV and was discovered due to a different relationship between RV and depending on the cycle phase, produced by the combination of geometrical effects (due to the butterfly diagram and the inclination) and the activity level variability along the cycle.

The existence of the hysteresis is a major result, because it comes as a limitation to the RV correction. We have shown that this hysteresis pattern is present in RV solar time series, as well as in other stars. Also, this result will help to better model the long-term RV due to activity.

We propose a new method which significantly improves the correction for long timescales (fraction of the cycle, i.e. suitable for planets in the habitable zone). The method is based on the physical description of the processes at the origin of the hysteresis (butterfly diagram related to dynamo processes and projection effects), which allows us to better constrain the relationship between the observables (RV and ) and therefore to better control the residuals, that is, to avoid adding spurious effects due to unphysical contributions. We note that Haywood et al. (2016) found that the unsigned magnetic flux (which is an observable that is currently very difficult to measure on solar-like stars other than the Sun) may be better correlated with long-term RV variation than : the difference may be due to a similar effect, with the unsigned flux corresponding to different projection effects. The new method proposed in the present paper needs to be further tested on all simulations to derive realistic detection rates for planets in the habitable zone around F6–K4 stars. After improvement and validation, these functions can then be made publicly available.

In addition to these new perspectives to detect low-mass planets in the habitable zone, we will also investigate the possibility of using this method to derive information on the spatio-temporal distributions of stellar activity, in particular in the direction (poleward or equatorward) of the dynamo wave, which would be extremely interesting since this has been very poorly constrained from observations so far. This should be feasible especially if the stellar inclination is well constrained.

Finally, the correction of the RV signal due to activity will be critical for the TESS and PLATO follow-ups. When these follow-ups aim at confirming planets observed by transit, this will concern stars seen close to edge-on statistically (unless they have a highly inclined orbit with respect to the equatorial plane of the star) and therefore with a significant hysteresis effect. For this reason, it is necessary to improve correction techniques taking this effect into account.

Acknowledgements

This work has been funded by the ANR GIPSE ANR-14-CE33-0018. We are very grateful to Charlotte Norris who has provided us the plage contrasts used in this work prior the publication of her thesis. The Ca K index was provided by the Sacramento Peak Observatory of the US Air Force Phillips Laboratory (https://www.nso.edu/uncategorized/ca-ii-k-line-monitoring-program/). This work was supported by the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES. This work was supported by the Programme National de Planétologie (PNP) of CNRS/INSU, co-funded by CNES. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Appendix A: RV–log R’_HK correlations

This appendix provides further details of the RV– correlations for our simulations: the correlation is computed either on whole time series (global correlations) or on subsets of data (local correlations), as discussed in Sect. 3. We then quantify these correlations for the Sun.

	Fig. A.1 Correlation between RV (due to activity only, SERIE 1, left, and including rvogs and rvinst, SERIE 2, right) and vs. average for different types of stars: F stars (upper panels), G stars (middle panels), and K stars (lower panels). The colour code is similar to that of Fig. 2. Only one point out of five is shown for clarity.
Open with DEXTER

	Fig. A.2 Upper panel: correlation between RV (due to activity only) and vs. RV jitterdue to convection only (rvconv) divided by total RV jitter (due to activity only). The colour code is similar to Fig. 2. Lower panel: same including rvogs and rvinst (for both the total RV jitter and the RV used to compute the correlation). Only one point out of five is shown for clarity.
Open with DEXTER

A.1 Impact of the parameters on the global correlation

The left panels of Fig. A.1 show the global correlation for RV due to activity only (rvspot2, corresponding to ΔT_spot2 as defined in Sect. 2), computed over each time series, as a function of the average activity level. The plots for ΔT_spot1 are very similar (not shown here). The correlation is usually high (above 0.7), but there is a drop in correlation for active F stars and high inclinations. This is due to the fact that for these stars the spot and plage contribution at the rotation period (which is not correlated with ) increases relative to the convection inhibition contribution: this decreases the correlation.

After adding the OGS signal as well as instrumental white noise (right panels), we observe a decrease of the correlation values for both quiet and active stars of all spectral types. The effect is particularly strong for K stars, despite the fact that the OGS signal is also decreasing when going to lower mass stars (see Paper I).

A.2 Relationship between correlation and convective RV component

Figure A.2 shows the correlation between RV and vs. the ratio between the RV jitter due to convection alone (which is the component which should be correlated to ) and the total RV jitter. With no-noise added (i.e. no OGS nor instrumental noise), there is a relatively good correlation between the two (0.92), but the relation is not entirely linear. The contribution of rvconv can be as low as 40–50% over our range of parameters, and the correlation between RV and , although always positive, can reach values of 0.2 (and even almost 0 in the presence of noise).

	Fig. A.3 Average short-term correlation between RV (including activity, rvogs and rvinst) and vs. global correlation. The colour code is similar to Fig. 2. Only one point out of five is shown for clarity.
Open with DEXTER

In the presence of noise (lower panel), the relationship is more complex. The global correlation remains good for two thirds of the simulations (66% still have a correlation above 0.7, while this was 95% for the no-noise case), but both ratios and correlations reach lower values (down to 0). Moreover, we see two regimes: one with a very good correlation and a linear relationship, and one below, where the correlation is degraded. This lower regime corresponds to F stars (37% of the F stars are in the second regime) with high inclinations (i.e. close to edge-on) compared to average. This is due to the fact that the correlation seems to be more strongly affected by inclination than the RV jitter ratio for these stars, leading to the observed shift toward lower correlations. This is probably due to the fact that in some cases, the addition of the OGS signal significantly degrades the correlation even though its contribution to the RV jitter is not major.

A.3 Short-term correlation

We have sofar considered the global correlation between RV and . This global correlation includes contributions from both short and long timescales. We now compute a short-term correlation, which is sensitive to the short (P_rot) timescales but not to the longer timescales: it is impacted by the convective blueshift inhibition in plages (correlated with ) and by the spot+plage signal, as well as the OGS signal when considered (not correlated with ). We compute the short-term correlation over time spans of 90 days: this is repeated for N consecutive intervals over the whole time series, with N between 35 and 59 depending on the length of time series. The average of the N values over each time series is then computed. The resulting averaged short-term correlation is shown in Fig. A.3.

The short-term correlation is related to the global correlation, however there is a large number of simulations with short-term correlations much lower than the global one. We also note that if the inclination of a given star is well known, it is possible to use this plot to relate the global and short-term correlations.

For RV including the OGS signal and the instrumental noise, the short-term correlations can vary within a range of 0.4–1 depending on the time series, as shown in Fig. A.4. For example, for an average short-term correlation of 0.6, the minimum varies between −0.2 and 0.4, andthe maximum between 0.75 and 0.95. For an average short-term correlation, the minimum is typically in the range −0.4 to −0.1 (and therefore negative) and the maximum between 0.45 and 0.8.

Fig. A.4 Upper panel: minimum value of short-term correlation (including activity, rvogs and rvinst). The minimum is computed over the N values of a given time series (see text). The colour code is similar to Fig. 2. Upper panel: same for the maximum of the short-term correlations. Lower panel: same for the maximum minus minimum of the short-term correlation. Only one point out of five is shown for clarity.

Open with DEXTER

A.4 The solar case

Figure A.5 shows the local correlation in the solar case for the two time series obtained by Meunier et al. (2010a,b) covering a complete cycle. They are described in Sect. 5.3. The average local correlation is similar (0.57–0.58) while the global correlation is much higher (0.91), which agrees with our simulations (see Figs. A.1 and A.3). The dispersion of the correlation values is respectively 0.20 and 0.27 and the local correlations cover a wide range of values as in our simulations.

	Fig. A.5 Local correlation between RV and vs. time (Julian days - 2 450 000) for RV reconstructed from observed solar structures (Meunier et al. 2010a) in black (stars) and from RV estimated from MDI/SOHO Dopplergrams (Meunier et al. 2010b) in red (diamonds). The horizontal lines correspond to the global correlations.
Open with DEXTER

Similar local correlations can be computed from the HARPS solar observation published by Milbourne et al. (2019), and are shown in Fig. A.6; the average of around 0.41 is slightly lower than for the previous time series. The correlation with the HMI/SDO RV that these latter authors reconstructed (no noise, spot, plage and convective blueshift inhibition only) is larger: on average around 0.63. However, if a noise of 1.2 m s⁻¹ is added for example, the average local correlation decreases (0.39, similar to the correlation with HARPS). The local correlations also show a wide dispersion, as in our previous results.

	Fig. A.6 Local correlation between RV and vs. time (Julian days – 2 457 222.5) for RV observed with HARPS (black stars), and reconstruted from HMI (raw in red diamonds, with noise added in green triangles), estimated from time series published by Milbourne et al. (2019). The horizontal lines correspond to the global correlations.
Open with DEXTER

References

All Tables

All Figures

	Fig. 1 Gain in RV jitter (when correcting for alone) vs. global correlation between RV and (including rvogs and rvinst).
Open with DEXTER
In the text

	Fig. 2 Gain in RV jitter in various conditions. Upper panel: gain vs. RV jitter (for simulations including rvogs and rvinst, and for ΔTspot2), for a correction with alone. The straight black line has a slope of 1, indicating points where the correction allows a residual RV jitter of 1 m s−1 to be reached. The colour code corresponds to inclination, from pole-on (i = 0°, yellow) to edge-on (i = 90°, blue), with light and dark orange corresponding to 20° and 30°, light and dark red to 40° and 50°, brown to 60°, and light and dark green to 70° and 80°. Middle panel: same for a correction with and a second-degree polynomial in time. Lower panel: gain with and time correction vs. gain with correction only. Only one point out of five is shown for clarity.
Open with DEXTER
In the text

	Fig. 3 Upper panels: amplitude of the polynomial in time vs. B−V (left) and (right). The black dots represent each individual time series, while the red points correspond to a binning in B− V and respectively. Lower panels: percentage of polynomial amplitudes higher than 0.5 m s−1 (stars), 1 m s−1 (diamonds), and 2 m s−1 (triangles) vs. B−V and respectively.
Open with DEXTER
In the text

	Fig. 4 Example of smoothed time series for a moderately active G2 star seen edge-on (left) and pole-on (right). First line: RV vs. time (solid line), and RV fitted from the variations (dashed line). Second line: vs. time. Third line: RV vs. , during the ascending phase of the cycle (red) and during the descending phase (green).
Open with DEXTER
In the text

	Fig. 5 Upper panels: hysteresis amplitude Chyst vs. B− V (left), vs. (middle), and vs. RV jitter (right), for F stars. The colour code is similar to that in Fig. 2. Middle panels: same for G stars. Lower panels: same for K stars. Only one point out of five is shown for clarity.
Open with DEXTER
In the text

	Fig. 6 Average hysteresis vs. RV jitter for three values of θmax: solar (black), solar +10° (red), and solar +20° (green). Each plot corresponds to a different inclination, from pole-on (0°, upper left panel) to edge-on (90°, lower right panel). All plots are on the same scale in m s−1.
Open with DEXTER
In the text

	Fig. 7 Relative velocity (solid line), RV deduced from the linear relation with Sacramento Peak S-index (dashed line) vs. time (left panels), and RV vs. S-index (right panels, in red for ascending phase and in green for descending phase). Upper panels: RVs reconstructed from observed structures (Meunier et al. 2010a). Lower panels: RVs reconstructed from MDI/SOHODopplergrams (Meunier et al. 2010b).
Open with DEXTER
In the text

	Fig. 8 Relative velocity vs. time (stars linked by solid black lines), and RV from linear dependence on (diamonds linked by orange dashed line) for each star. The second plot for each star shows the hysteresis pattern, with the ascending phase in red and descending phase in green.
Open with DEXTER
In the text

	Fig. 9 Average vs. cycle phase for pole-on (solid line) and edge-on (dashed line) configurations for simulation shown in Fig. 4.
Open with DEXTER
In the text

	Fig. 10 Relative velocity vs. for simulation similar to Fig. 4, seen pole-on (upper panel) and edge-on (lower panel), but with constant latitude over time.
Open with DEXTER
In the text

	Fig. 11 First panel: example of correcting function from catalogue (see Sect. 6.4) vs. phase for edge-on configuration (left) and pole-on configuration (right): chromospheric emission (solid line), RV (dashed line), and ratio between RV and chromospheric (thick dotted-dashed line), all normalized to 1 for clarity. Second panel: RV vs. chromospheric emission for the simulation shown in Fig. 4, colour-coded according to the cycle phase: red for the ascending phase, green for the descending phase. The orange and brown curves correspond to the same phases respectively but for the model derived from the catalogue functions. The black dashed line is the linear fit between RV and . Third panel: RV residual after standard correction (black dashed line) and after the new correction (orange and brown solid lines as in the second panel). Fourth panel: RV residual vs. time after standard correction (dashed line) and after the new correction (solid line).
Open with DEXTER
In the text

	Fig. 12 Distribution of different rms ratios, for inclinations of 0° (orange), 30° (red), 60° (brown), and 90° (green). Panel A: rms0/rmsstand. Panel B: rms0/rmsnew. Panel C: rmsstand/rmsnew. Panel D: rmsnew/rmsopt. Panel E: rmsnew/rmsext. Panel F: rmsext/rmsopt.
Open with DEXTER
In the text

	Fig. A.1 Correlation between RV (due to activity only, SERIE 1, left, and including rvogs and rvinst, SERIE 2, right) and vs. average for different types of stars: F stars (upper panels), G stars (middle panels), and K stars (lower panels). The colour code is similar to that of Fig. 2. Only one point out of five is shown for clarity.
Open with DEXTER
In the text

	Fig. A.2 Upper panel: correlation between RV (due to activity only) and vs. RV jitterdue to convection only (rvconv) divided by total RV jitter (due to activity only). The colour code is similar to Fig. 2. Lower panel: same including rvogs and rvinst (for both the total RV jitter and the RV used to compute the correlation). Only one point out of five is shown for clarity.
Open with DEXTER
In the text

	Fig. A.3 Average short-term correlation between RV (including activity, rvogs and rvinst) and vs. global correlation. The colour code is similar to Fig. 2. Only one point out of five is shown for clarity.
Open with DEXTER
In the text

	Fig. A.4 Upper panel: minimum value of short-term correlation (including activity, rvogs and rvinst). The minimum is computed over the N values of a given time series (see text). The colour code is similar to Fig. 2. Upper panel: same for the maximum of the short-term correlations. Lower panel: same for the maximum minus minimum of the short-term correlation. Only one point out of five is shown for clarity.
Open with DEXTER
In the text

	Fig. A.5 Local correlation between RV and vs. time (Julian days - 2 450 000) for RV reconstructed from observed solar structures (Meunier et al. 2010a) in black (stars) and from RV estimated from MDI/SOHO Dopplergrams (Meunier et al. 2010b) in red (diamonds). The horizontal lines correspond to the global correlations.
Open with DEXTER
In the text

	Fig. A.6 Local correlation between RV and vs. time (Julian days – 2 457 222.5) for RV observed with HARPS (black stars), and reconstruted from HMI (raw in red diamonds, with noise added in green triangles), estimated from time series published by Milbourne et al. (2019). The horizontal lines correspond to the global correlations.
Open with DEXTER
In the text