Activity time series of old stars from late F to early K. IV. Limits of the correction of radial velocities using chromospheric emission

Inhibition of the convective blueshift in active regions is a major contribution to the radial velocity variations, at least for solar-like stars. A common technique to correct for this component is to model the RV as a linear function of chromospheric emission, because both are strongly correlated with the coverage by plages. This correction is not perfect: the aim of the present study is to understand the limits of this correction and to improve it. We investigate these questions by analysing a large set of synthetic time series corresponding to old main sequence F6-K4 stars modelled using a consistent set of parameters. We focus here on the analysis of the correlation between time series, in particular between RV and chromospheric emission on different timescales. We also study the temporal variation for each time series. Inclination strongly impacts these correlations, as well as additional signals (granulation and supergranulation). Although RV and LogR'HK are often well correlated, a combination of geometrical effects (butterfly diagrams related to dynamo processes and inclination) and activity level variations over time create an hysteresis pattern during the cycle, which produces a departure from an excellent correlation: for a given activity level, the RV is higher or lower during the ascending phase compared to the descending phase of the cycle depending on inclination, with a reversal for inclinations about 60 deg from pole-on. This hysteresis is also observed for the Sun and other stars. This property is due to the spatio-temporal distribution of the activity pattern and to the difference in projection effects of the RV and chromospheric emission. These results allow us to propose a new method which significantly improves the correction for long timescales, and could be crucial to improving detection rates of planets in the habitable zone around F6-K4 stars.


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 lowmass 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 log R ′ HK , 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 Send offprint requests to: N. Meunier 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 log R ′ HK 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 log R ′ HK and this polynomial.
(2) The correction using the log R ′ HK 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 to 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 A&A proofs: manuscript no. ms35348fin I, which also shows how the log R ′ HK , representing the 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-R' HK 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 log R ′ HK time series across our grid of parameters. Subsequently, we present a study of the typical gain which is obtained when performing standard correction using log R ′ HK and a polynomial in time. In Sect. 5, we identify a limiting element in the log R ′ HK 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.

Model and parameters
Our model, described in detail in Borgniet et al. (2015) for the Sun and in Paper I for other stars (F6 to 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 1 ), 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 hour, and an instrumental white noise represented by a Gaussian noise with an amplitude of 0.6 m/s (rvogs). Here, log R ′ HK 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 2019c) 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 hours). The time series have a maximum length of 15 years, 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.

HK
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-R' HK 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 2 between RV and log R ′ HK , because it can easily be derived from observational data and gives some clues on what to expect from the log R ′ HK correction. Typical values are shown in detail in Appendix A.1. The correlation between RV and log R ′ HK 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 log R ′ HK 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 log R ′ HK includes contributions from both short and long timescales. The shortterm correlations presented in Appendix A.3 are 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, we observe 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 shortterm 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

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=α+βlog R ′ HK . 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=α+βlog R ′ HK +γt+δt 2 (Dumusque et al. 2017).

Gain while correcting from log R ′ HK and polynomial in time
It is interesting to see whether it is easier to correct for activity using the log R ′ HK -RV correlation when the correlation between RV and log R ′ HK is high. Figure 1 shows the gain defined at the beginning of Sect. 4 versus the log R ′ HK -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 log R ′ HK 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 RV jitter (for simulations including rvogs and rvinst, and for ∆Tspot2), for a correction with log R ′ HK alone. The straight black line has a slope of 1, indicating points where the correction allows a residual RV jitter of 1 m/s 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 log R ′ HK and a second-degree polynomial in time. Lower panel: Gain with log R ′ HK and time correction vs. gain with log R ′ HK correction only. Only one point out of five is shown for clarity.
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 log R ′ HK and a second degree polyno-Article number, page 3 of 15 A&A proofs: manuscript no. ms35348fin mial in time (which is often done to 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 log R ′ HK and polynomial in time correction vs. gain with log R ′ HK 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 log R ′ HK , we use a model which is not physical.

Discussion on the polynomial function fitted to the RV time series
Given the difference in gain observed for certain simulations when correcting for log R ′ HK 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 log R ′ HK 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 log R ′ HK (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 log R ′ HK 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.

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.

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 ∆Tspot 2 , and include activity, OGS (smoothed over 1 hour), 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 log R ′ HK were used. The lower panel shows RV versus log R ′ HK : 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 log R ′ HK 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 log R ′ HK , 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.

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 log R ′ HK levels (corresponding to the range covered by each time series). The RV is averaged over each log R ′ HK 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, log R ′ HK , 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 (edgeon) 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 log R ′ HK . This is detailed and discussed in Sect. 6, where a new correction method is proposed.

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.

The solar case
Our simulations provide evidence for the presence of an hysteresis pattern between the long-term RV and log R ′ HK 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 in 2019). The reconstruction of the A&A proofs: manuscript no. ms35348fin

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 Notes. The spectral type is from the CDS (https://simbad.ustrasbg.fr/simbad/). "E" indicates an hysteresis pattern similar to our edge-on configurations, "P" similar to our pole-on configurations, "M" a mixed pattern, and "O" no hysteresis.
Article number, page 6 of 15 Meunier et al.: Activity time series of old stars from late F to early K. IV.  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.

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 log R ′ HK , and illustrate its performance on a subset of G2 star simulations.

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 Article number, page 7 of 15 A&A proofs: manuscript no. ms35348fin the velocity field (µ), a contrast dependence on µ, and the centerto-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 a small range of values: for a given activity level, the RV signal is lower during the ascending phase of the cycle. This be- haviour 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 strong effect 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 relation-A&A proofs: manuscript no. ms35348fin ship between RV and log R ′ HK , although showing little dispersion, is not strictly linear and some dispersion is still present.

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 log R ′ HK 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 20930 functions. Figure 11 shows an example of such a function (upper panels).

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 poleon 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 ∆Tspot 1 . We compare the rms before correction (rms 0 ), 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.

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 nonlinearity observed in Fig. 10. After correction of this curvature using a second-degree polynomial in log R ′ HK 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 log R ′ HK 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 Notes. Median of the gains between the different correction levels for four inclinations between pole-on (0 • ) and edge-on (90 • ). These correspond to the distributions shown in Fig. 12.
RV=α(1+βlog R ′ HK +γ(log R ′ HK ) 2 )(1+δt+ǫt 2 ) where the polynomial in log R ′ HK 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).

Conclusion
We analysed a large number of simulated time series of stellar activity covering the spectral types F6 to 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 log R ′ HK . The detailed study of the correlation between RV and log R ′ HK , and the gains in RV jitter which can be obtained using the typical log R ′ HK 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 log R ′ HK limited (as shown in previous papers Meunier & Lagrange 2013;Dumusque et al. 2017;Meunier & Lagrange 2019a), 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 log R ′ HK was discovered due to a different relationship between RV and log R ′ HK 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 log R ′ HK ) 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 log R ′ HK : 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 spatiotemporal distributions of stellar activity, in particular in the direction (poleward or equatorward) of the dynamo wave, which A&A proofs: manuscript no. ms35348fin  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. Correlation between RV (due to activity only, SERIE 1, left, and including rvogs and rvinst, SERIE 2, right) and log R ′ HK vs. average log R ′ HK 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. Appendix A.2: Relationship between correlation and convective RV component Figure A.2 shows the correlation between RV and log R ′ HK vs. the ratio between the RV jitter due to convection alone (which is the component which should be correlated to log R ′ HK ) 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 log R ′ HK , although always Correlation between RV (due to activity only) and log R ′ HK vs. RV jitter due 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. positive, can reach values of 0.2 (and even almost 0 in the presence of noise).
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.

Appendix A.3: Short-term correlation
We have so far considered the global correlation between RV and log R ′ HK . 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 log R ′ HK ) and by the spot+plage signal, as well as the OGS signal when considered (not correlated with log R ′ HK ). 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 to 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, and the 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.  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. Local correlation between RV and log R ′ HK vs. time (Julian days -2457222.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.