The effects of granulation and supergranulation on Earth-mass planet detectability in the habitable zone around F6-K4 stars.

We first define the solar values we consider in this study. The time series are derived from power spectra following Harvey (1984) for granulation and supergranulation and following the shape of the envelope of the oscillations from Kallinger et al. (2014), as in Papers I and II. This method has the advantage of allowing us to produce a large amount of very long time series. We showed in Meunier et al. (2015) that the shape proposed by Harvey (1984) was well adapted, even at low frequencies: therefore, we use the parameters found in Meunier et al. (2015). The choice of a one-hour binning is similar to what we chose in Paper I and corresponds to the timescales where the RV jitter due to granulation reaches an inflexion point Meunier et al. (2015): binning over a longer duration is not efficient enough to reduce this jitter further and so, this binning time is used to filter granulation out best.

For granulation, in the majority of our study, we use an rms (root-mean-square) of 0.83 m/s before averaging (i.e. 0.39 m/s after averaging over one hour), hereafter GRAhigh, which stands as our reference value, provided by our simulations of about 1 million granules on the disk at any given time in Meunier et al. (2015). As discussed in Paper I, such simulations were based on realistic properties of granules (derived from hydrodynamical simulations of Rieutord et al., 2002), which are known to reproduce realistic line profiles Asplund et al. (2000). However, lower values were derived from the observation of two specific spectral lines: about 0.32 m/s by Elsworth et al. (1994) from the Potassium line at 770 nm and 0.46 m/s from the Sodium doublet at 589 nm by Pallé et al. (1999). More recently, the residuals on timescales lower than one day obtained by Collier Cameron et al. (2019) on solar integrated RV times series and covering the whole spectrum obtained by HARPS-N are also of the order of 0.40 m/s (when averaging over typically five minutes). Similar amplitudes have been obtained by Sulis et al. (2020) using MHD simulations. The difference between these estimates and the results of Meunier et al. (2015) may be due to some subtle effects in the centre-to-limb dependence which are not taken into account in Meunier et al. (2015), but also to the fact that observations were made in a single lines which may not be representative of the whole spectrum. For that reason, a twice lower level (hereafter GRAlow) with respect to our reference level will also be considered in mass characterisations and blind tests in Sect. 4.3 and 4.4. We note that Cegla et al. (2019) obtained very low rms RV for granulation using a reconstruction based on MHD simulations of the solar surface, around 0.1 m/s. The reason for this discrepancy is not clear at this stage, although it may be due the fact that strong vertical magnetic fields were used.

Concerning supergranulation, Meunier et al. (2015) provide a large range of possible values based on our current knowledge of these flows. Here, we consider two values, their median level (0.7 m/s, hereafter SGmed), and their lower level (0.27 m/s, hereafter SGlow), as in Meunier & Lagrange (2019b): these are in agreement with typical amplitudes obtained for a few stars by Dumusque et al. (2011). The median level is also close to the rms found by Pallé et al. (1999) for the Sun, with 0.78 m/s for the Sodium doublet lines. Because of the longer timescales of supergranulation, the rms RV is almost the same after the 1 hour averaging. The amplitude of the oscillations is derived from Davies et al. (2014), as in Paper I. The time scale is the same one obtained in Meunier et al. (2015) as in Paper I, that is, 1.1 10⁶ s.

We mainly use five types of time series throughout the paper: high level of granulation alone (GRAhigh), supergranulation alone (SGmed, median level, and SGlow, low level), all contributions for oscillations, a high level of granulation, and median supergranulation (ALL_{GRAhigh,SGmed}) or low supergranulation (ALL_{GRAhigh,SGlow}). In the following, ALL always represents the superposition of oscillations, granulation, and supergranulation. The other three configurations (GRAlow alone, ALL_GRAlow,SGmed, ALL_GRAlow,SGlow) are mostly be considered for the mass characterisation and blind tests to provide a complete view. The configuration ALL_{GRAhigh,SGmed} was used in combination with magnetic activity in Paper II. The contribution attributed to any of these combinations is referred to as the OGS (for oscillations, granulation, supergranulation) signal in the following. The oscillations are not studied alone here because we consider one-hour averages and they are well averaged out Chaplin et al. (2019) at such timescales: they did not prevent us from obtaining excellent detection rates when considered independently (Paper I).

We considered seven spectral types covering the F6-K4 range, that is, F6, F9, G2, G5, G8, K1, and K4. The amplitudes of the different components were scaled to G2 stars (i.e. solar values from the previous section) as in Paper II (previous section). We recall them here in brief. Granulation parameters are scaled from G2 stars to other spectral types using results from Beeck et al. (2013). Oscillation parameters are scaled using laws from Kjeldsen & Bedding (1995), Samadi et al. (2007), Bedding & Kjeldsen (2003), Kippenhahn & Weigert (1990), and Belkacem et al. (2013)¹¹1There are other scaling laws, for example by Yu et al. (2018), but they are not very different. Because the oscillations are strongly averaged in this work, this choice is not critical.. Supergranulation is scaled following the granulation scaling, assuming supergranulation is strongly related to granulation properties Rieutord et al. (2000); Roudier et al. (2016), including the time scale, which can differ by up to about 20%, so that the impact should be small.

All time series were produced for a ten-year period of duration with a 30-second time step and are then binned over one hour. We then selected one such point per night. Examples of time series (subsets over short periods) are shown in Appendix A for F6, G2, and K4 stars, as well as examples of the power functions versus frequency. In addition to this full sample of 3650 nights, we consider several other configurations with a gap of four months per year to simulate the fact that a star can usually not be observed all year long. Then $N_{\rm obs}$ nights were randomly selected out of the remaining nights over the ten-year duration. Each realisation of this selection for a given value of $N_{\rm obs}$ corresponds to a different sampling. We use $N_{\rm obs}$ = 180, 542, 904, 1266, 1628, 1990, and 2352 nights with the four-month gap each year), and 3650 nights (no gap), leading to a total of eight configurations. In Paper I, we found that using a random selection or considering packs of adjacent nights did not lead to significant differences. In addition, Burt et al. (2018) tested different ways of building the sampling for magnetic activity time series and found that the random sampling was optimal (the uniform sampling was not extremely different however). Testing of additional sampling configurations is presented in Sect. 5.

Figure 1 summarises the rms RV versus spectral type for the eight configurations of OGS time series used in this paper. We note a general decrease towards lower mass stars. When considering all components and spectral types, the RV jitter varies between 0.28 and 0.9 m/s typically when considering GRAhigh. The dashed lines show the levels when the granulation level is divided by 2 (GRAlow). In this case the granulation rms varies between 0.22 and 0.1 m/s, and when combined with the low level of supergranulation it varies between 0.2 and 0.37 m/s. We note that even for such a large number of points, there is little dispersion in RV jitter from one realisation to the next.

We focus our analysis on low-mass planets orbiting in the habitable zone of their host stars. We define the limits of the habitable zone as a function of spectral type as in Meunier & Lagrange (2019a), following Kasting et al. (1993), Jones et al. (2006), and Zaninetti (2008). We consider three typical orbital periods, corresponding to the inner side (PHZ_in), the middle (PHZ_med), and the outer side (PHZ_out) of the habitable zone: the resulting orbital periods vary between 409-1174 days for F6 stars to 179-501 days for K4 stars. Furthermore, we consider only circular orbits, for simplicity.

Most of the computations are carried out with projected masses of 1 and 2 M_Earth. For inclinations higher than 40-50^∘, the performance obtained with these masses is representative of this whole range of inclinations, while for lower inclinations the performance should be significantly worse than the one presented in this paper. Therefore, additional blind tests, presented in Sect. 4.4, are also performed, considering a distribution of inclinations between 0^∘ and 90^∘ when building the data set and with the assumption that the orbital plane is the same as the stellar equatorial plane. In the case of a transit follow-up using RV to characterise the mass, however, the projected mass can be considered to be the true mass.

In subsequent sections, the analysis of the time series is made using two complementary approaches (i.e. two test statistics), which are then compared. The steps are as follows: (i) we analyse the periodograms²²2We use the Lomb Scargle periodogram with no normalisation to be able to compare powers between different types of contributions. They are computed between 2 and 2000 days. of the time series, and compute the maximum amplitude around the considered PHZ (frequential analysis, computed in 0.9–1.1 PHZ range); (ii) we fit the planetary signal, considering a period guess corresponding to the period of this peak with maximum amplitude (temporal analysis) or of interest (PHZ) depending on the case. This fit is made using a $\chi^{2}$ minimization.

We first consider the detectability of such exoplanets in the presence of the stellar contribution defined in the previous section. Because we consider synthetic time series, we can study them with the certainty that there is no planet present in the signal. As a consequence, we can estimate a true level of false positive (FP) for a given test statistics (frequential and temporal analysis) and for a given probability (e.g. 1%), and it is then possible to compute detection rates for a given planet (on the time series where the planet have been added), considering this level of false positives. The method we apply to determine the FP is described in Sect. 2.3. Once we have determined a true FP level corresponding to a certain percentage of false positives, and a detection rate for a given mass, we can also determine which mass corresponds to a good detection rate (e.g. 95%), which provides a detection limit. This approach is explored in Sect. 3.

From the point of view of the observer, however, the determination of the true level of false positive due to a given stellar contribution is not possible because it is not possible to know if it includes other, additional signals (of a planet for instance) and because we have only one realisation of the signal. This is why the analysis of observed time series always relies on other methods, such as the use of false alarm probability levels using bootstrap analysis, although this approach makes assumptions on the signal which may not be correct, as pointed out by Sulis et al. (2017a). In a second step, we therefore test this type of approach and compare it with the one based on the true false positive level. We also compare the detection limits based on the periodogram analysis proposed by Meunier et al. (2012), the local power analysis (LPA) method with the true detection limits. A blind test is implemented to estimate the detection rates and false positive levels and compare them with the true ones. This approach is explored in Sect. 4.

The latter issue, also studied in Sect. 4, concerns the performance with regard to mass characterisations of planets detected by transit in photometric light curves. In this case, we consider that the planet presence is confirmed,meaning that the transits do not require any validation using RV observations. There is, therefore, no issue with false positives in this case and we also know its orbital period and phase with very good precision from the transit. We can then fit the RV amplitude due to the planet (temporal analysis) at this orbital period to determine the precision for the mass characterisation.

For each OGS signal and spectral type, statistically half of the realisations of the time series remain unchanged while a planet is added to the other half. The analysis of each time series allows to determine whether a planet is detected on not. In a second independent step, we determine the level of false positives and the detection rate for each set of simulations, by comparing the outputs with what was actually injected or not. We focus our analysis on one of the $N_{\rm obs}$ values (1266 points), which corresponds to good conditions, but with still a reasonable rate of observations in future dense monitorings.

The fitting challenge implemented in Dumusque et al. (2017) which focusses on stellar magnetic activity defined several detectability criteria. We use similar criteria and terminology with a few modifications: 1) We decide whether there is a detection or not using a binary choice, but since there is no further comparison with activity indicator for example, there is no intermediate case; 2) False positives are counted separately for realisations with an injected planet and with no planet; 3) The identification of the planet in Dumusque et al. (2017) was based on whether the amplitude K and the period P corresponded to the injected planet, while here we use only the period as a criterion, because given the dispersion in mass this criteria would be quite subjective and can be used in a second step if necessary. The different categories are summarised in Table 2. False positives and missed or wrong planets can bias statistics on exoplanets: their effect is also indicated in Table 2. We also note that the detection criteria in Dumusque et al. (2017) was not the same in all methods as it depended on the fitting method, and may be different from ours.

The first step of the procedure consists of building the data sets. For each configuration (one spectral type, and one of the eight OGS configurations), we consider 200 realisations of the OGS signal and sampling. Computations are made for 1266 points and 1 M_Earth unless otherwise noted (figures for other values and approaches such as including a distribution in inclinations are shown in Appendix C). Based on a random variable, on average, half of the realisations remain unchanged, while a planet is added to the other half. The planet has the following properties: the orbital period P_orb is chosen randomly in the PHZ_in-PHZ_out range (i.e. we consider the habitable zone globally, using a uniform distribution), and the phase is chosen randomly in the [0-2$\pi$] interval. The projected mass is equal to 1 M_Earth (projected mass, see discussion in Sect. 2.2.1) unless noted otherwise: these blind tests serve as our reference. Figures corresponding to other masses are shown in Appendix C, along with blind tests that include inclination distribution.

For each configuration (one spectral type, and one of the 8 OGS configurations), each of the 200 realisations of the time series are analysed as follows. The FAP at the 1% level is computed (using 200 realisations of the bootstrap, which we checked does not give very different results from a larger number of realisations). The periodogram of the time series is computed and the highest peak is identified (in the range 2-2000 d). If the amplitude of the peak is lower than the FAP, then we establish that there is no detection, whereas if the amplitude of the peak is higher, we consider this to be a detection. In this latter case, a sinusoid (we recall that we consider only circular orbits here) is fitted, with the period fixed to the peak period, to obtain the mass.

We note that the conditions are different from the theoretical results presented in Sect. 3. In Sect. 3, each computation was focusing on the behaviour at a given period (for example the middle of the habitable zone) and, therefore, on the power at this particular period or the mass corresponding to a fit at this period. Here, we address a different question, since we do not focus on a given period: we place ourselves at the point of view of an observer and we do not know where the planet is injected, that is, we consider the whole 2-2000 day range and not only the habitable zone. The analysis can even lead to a (wrong) detection outside of the period range where the planet is injected. This can then induce a higher rate of false positives (unless the criteria to make the detection is much higher than the true false positive level).

In a second independent step, we compare the results with the input parameters: this allows us to associate one of the categories of Table 2 to each realisation. The decision algorithm is shown in Fig. 15. To define whether a peak is attributable to the true planet or not, we use the criterion $|P_{\rm peak}-P_{\rm true}|<0.1P_{\rm true}$ to determine if the planet is the correct one (see next section).

The outputs of each blind test are mainly the properties of the fitted planet parameters when detected and the percentages corresponding to the different categories defined in Table 2. We first focus on the period since a criterion on the period obtained during the analysis must be defined to assign each realisation to one of the categories. Figure 16 shows the distribution of the difference between the periods provided by the analysis and the true periods over all realisations (i.e. all OGS configurations and all realisations with an injected planet), independently of the significance of the peak. The realisations outside this range correspond mostly to peaks found at low periods, with a maximum of the distribution in the 20-30 days range, as shown in the lower panel of Fig. 16: 95% of those peaks are at periods below the true orbital period and many of them are, in effect, below the FAP. In practice, the width of the peak varies with the period, and a threshold of 10% of the period allows us to separate the peak from outliers.

Table 3 shows a few examples of percentages, for G2 and K4 stars and a subset of OGS (GRA, ALL_{GRAhigh,SGmed}, ALL_{GRAhigh,SGlow}). Ideally, we would like to obtain 100% on the first two lines and 0% on the other lines. The categories correspond to Table 2, some of them being regrouped. For example, the percentage of bad planet detections (i.e. the global false positive rate) corresponds to planets detected when none was injected, added to the planet detected with a wrong period. For G2 stars and granulation, the recovery rate is very good when no planet is injected but lower when there is an injected planet: most of the lost planets correspond to peaks below the FAP. The performance is much poorer for ALL_{GRAhigh,SGmed}, with very low detection rates when a planet is injected and high false positive level. Even for ALL_{GRAhigh,SGlow}, the recovery rate when a planet is injected is only 35%. For K4 stars, performance is perfect of granulation and very good for ALL_{GRAhigh,SGlow}. For ALL_{GRAhigh,SGmed}, the detection rate is only 67%, however.

Figure 17 summarises the percentages for all configurations (1266 points, 1 M_Earth). The good recovery rates are shown on the left-hand side panels. When no planet is injected (black curves), they are very good for GRAhigh and GRAlow, and above 80-90% when added to SGlow. They are strongly degraded in other configurations, for all spectral types (and more so for high mass stars). The detection rates when a planet is injected (green curves) are good for all stars for GRAlow only, and for K stars and sometimes G stars for GRAhigh, SGlow, ALL_{GRAhigh,SGlow}, and ALL_GRAlow,SGlow (the threshold depends on the configuration) but strongly degraded for all other cases. It could seem surprising that the performance is better when considering ALL_{GRAhigh,SGlow} compared to SGlow alone (no injected planet): this is likely due to the fact that when adding the GRAhigh signal, even though the rms is increased, the power spectrum is then more similar to the GRAhigh power spectrum corresponding to good performance in the habitable zone.

The green dotted lines correspond to the detection rate obtained for the theoretical false positive level of 1% (Sect. 3), which is to be compared to the green solid line observed in the blind test. The two estimations are sometimes similar, corresponding to the FP that is very close to the FAP (Sect. 3), while in other cases, the blind test detection rates are lower than expected from the theoretical false positive level due to the difference between the FAP and the true FP. There is, therefore, a complex relationship between the theoretical results and the detection rates derived from the blind tests. We conclude that the FAP provides a detection rate which corresponds to a different false positive level from the one expected (i.e. in our case, diverging by 1%).

The right-hand panels in Fig. 17 show the bad recovery rates. When a planet is injected, the bad recoveries (dashed black line) naturally serve as the complement of the green curve from the left panels. It represents a wide variety of situations: it is sometimes dominated by the missed planet (bad period and below the FAP, in blue), sometimes by the rejection of true planet (in orange); globally, that is, because peaks are below the FAP and sometimes because the highest peak is above the FAP but does not correspond to the planet (in brown, i.e. a false positive). We note that the false positive rate when a planet is injected (brown) is different from the false positive when no planet is injected (in red, completary to the black curves on the left-hand side panels) for supergranulation (especially SGlow) alone, but it is similar when granulation and supergranulation are superposed.

For GRA, ALL_{GRAhigh,SGmed}, and ALL_{GRAhigh,SGlow}, the red and brown curves are similar, that is, the percentage of false positives is the same whether a planet (of 1 M_Earth) is injected or not. However, the situation is different for SGmed and SGlow because when no planet is injected, the percentage of false positives is the same, even though they have very different rms RV. This is due to the fact that here the comparison of the power is made with the FAP. Because the shape of the power spectrum is the same between SGmed and SGlow, and because the FAP values are scaled with the rms, both power and FAP increases from SGlow to SGmed in a similar manner and the percentage of false positives is similar. In the case of ALL_{GRAhigh,SGlow}, the signal is dominated by GRAhigh, hence, a level that is similar betwen GRAhigh and ALL_{GRAhigh,SGlow}, while the situation is intermediate for ALL_{GRAhigh,SGmed}. For GRAlow, the rate of false positives is very small in all cases. However, when added to supergranulation (either SGmed or SGlow), the latter dominates, and rates are very similar to those obtained when combining with GRAhigh, only slightly lower.

The level of false positives here may be large because our analysis is too simplistic. When a peak is detected above the FAP, we should test the robustness of the detection to determine whether the peak is stable or not for example. More sophisticated methods will have to applied in this area in the future (see Sect. 6).

Another representation of these results is shown in Fig. 18, showing the percentage of false positives (sum of the two contributions described above) versus the detection rate (computed on the cases with an injected planet), which is similar to a ROC curve (but where each point corresponds to a spectral type). Each curve corresponds to one of the OGS configurations. Ideally, we would like points to be in the lower right corner. Points at the top have a high false positive level and points on the left correspond to poor detection rates. If we compare the global level of false positive here and the rms for each type of OGS configuration, we see that there is not a direct correspondence, because a granulation-like signal provides better performance due to their more suitable power spectrum (for a given rms). High-mass stars are to the left of each curve and lead to high rates of false positives and low detection rates except for granulation alone (for GRA_low all points are in the lower right corner), and to a lesser extent ALL_{GRAhigh,SGlow}. We also note that the highest level of false positives is obtained for SG alone. However, when granulation is added to supergranulation, the rms increases, but the level of false positive decreases because the shape of the power spectrum is closer to the granulation shape, leading to better performance: This explains why the level of false positive is higher when SGmed and GRAlow are superimposed (dashed green curve) compared to SGmed and GRAhigh (solid green curve), that is, closer to the SG behavior (large false positive rates) even though the rms is lower.

Additional configurations are tested in Appendix C.1 (180 points only) and C.2 (2 M_Earth). The performance for 180 points is very poor. The level of false positives is quite low, which can be explained by the results shown in Sect. 4.1: here, the FAP overestimates the true false positive level and, therefore, there are few peaks above the FAP. The detection rates are very low, however. On the other hand, the performance is much better for a 2 M_Earth planet compared to a 1 M_Earth, although it is not perfect in all cases: for F and early G stars, the detection rates reach values below 50% when supergranulation is high.

We also implemented a similar blind test, but in which 1 M_Earth or 2 M_Earth are the true planet mass. We assume that the orbital plane is similar to the equatorial plane and take the distribution of stellar inclination into account. We expect slightly lower detection rates than before (for cases with injected planets), which is indeed observed as shown in Appendix C.3 and C.4. Figure 19 shows the average of the rates over all spectral types for each OGS configuration, without taking inclination into account (previous results) and, conversely, taking it into account. The detection rates are slightly lower when considering inclination (i.e. the true mass), typically by a different of about 12-13 points on the percentage. The difference is mostly due to the larger amount of missed planet when the mass is the projected mass only.

Finally, we compute the LPA detection limits (see Sect. 4.2 for the definition): with an injected planet with a mass of 1 M_Earth, we want the LPA detection limit (M_lpa) to be higher than 1 M_Earth. We compute ten values of M_lpa over the habitable zone, which are then averaged together for each spectral type. The average M_lpa and the percentage of realisations where M_lpa is higher than 1 M_Earth. In all cases, M_lpa is indeed above 1 M_Earth, and the percentage above 70%, which is in agreement with expectation. When no planet is injected on the other hand, we want M_lpa to be as low as possible. For SGmed and ALL_{GRAhigh,SGmed}, they are above 1 M_Earth for F6-G8 stars, so that in those cases, the exclusion of the presence of low mass planet (below 1 M_Earth) is not possible. This is strongly related to the performance in terms of detection rates described above. For all other configurations (OGS, spectral types), they are always below 1 M_Earth. We conclude that the LPA provides results that are consistent with the presence of the injected planet.

In this section, we compute the K/N criterion proposed in Dumusque et al. (2017) and defined as $K_{\rm pl}\sqrt{N_{\rm obs}}$/RV_rms, where $K_{\rm pl}$ is the amplitude of a planetary signal in RV (for a given mass, period, host star), $N_{\rm obs}$ is the number of observation, and RV_rms is the RV jitter⁷⁷7In the original definition in Dumusque et al. (2017), the RV jitter is computed after correction with a linear correlation with chromospheric emission index and trend in time, which is not done here because it is irrelevant to the type of signal we consider.. K/N is used by Dumusque et al. (2017) as a criterion for estimate the quality of recovery rates. Therefore, we compute this practical criterion for a 1 M_Earth planet and compare it to the detection rates obtained previously for the same planet mass (cases with injected planet). The results are shown in Fig. 20. We find a very clear relationship between the two: all OGS configurations and spectral types lie along the same curve with very little dispersion, so the criterion is adequate to describe the detection rate in these conditions. Detection rates better than 50% correspond to K/N above $\sim$7, and K/N must be above $\sim$9 to reach detection rates better than 95%. This is very similar to the rough threshold between bad recoveries and good recoveries of $\sim$7.5 in Dumusque et al. (2017), who focused on magnetic activity. On the other hand, there is not a one-to-one relationship between this criterion and the false positives, as the different OGS configurations correspond to different levels, as shown on the lower panel of Fig. 20.

Although the curve for a given number of points, $N_{\rm obs}$ , and the mass are well-defined, it is, in fact, very dependent on the configurations. For example, for a lower number of points (see Fig. 33 in Appendix C.1 for 180 points) the curve is very different: the curve is also well-defined, but for a similar K/N, the detection rates are lower than for 1266 points. The same is observed for 2 M_Earth, with the 50% level reached at lower K/N values compared to a 1 M_Earth planet. Thus, the criterion is not universal. We then consider the performance as a function of the number of points in Sect. 5.