© S. Sulis et al. 2020

1 Introduction

At the time of this writing, 880 extrasolar planets have been discovered so far by the radial velocity (RV) technique¹. In this sample, 52% of these consist of planets that are more massive than Jupiter and 14% that have a mass inferior to 10 Earth-masses (M_⊕). Since the detection of HD215152c (Mayor et al. 2011), only 19 planets have been found with a mass ≤2 M_⊕ and all of the latter are short-period (≤50 days) planets orbiting stars that are less massive than the Sun.

Indeed, detecting planets is easier around low-mass stars (as the ratio of the planet to stellar masses is higher) and, thus, a first strategy consists in monitoring cool M dwarfs to increase the detection probability. This has been the purpose of recent surveys with spectrographs such as CARMENES (Quirrenbach et al. 2014) and SPIRou (Donati et al. 2017). On the other hand, new instruments such as ESPRESSO (Pepe et al. 2010) and EXPRES (Jurgenson et al. 2016) have been developed to ensure the long-term stability that is needed to detect signals of terrestrial planets orbiting main sequence G-dwarf stars (with an amplitude around 10–30 cm s⁻¹).

However, detecting planet signatures at the cm s⁻¹ level remains challenging as spurious Doppler shifts of various origins may dominate the RV series. The activity at the surface of the host star is one of the main sources generating changes in depth, width, and asymmetries of the absorption lines. Disentangling the planetary signal from the stellar activity “noise” is an active research topic (see e.g., Desort et al. 2007; Aigrain et al. 2012; Haywood et al. 2014; Lagrange et al. 2010; Meunier et al. 2017a; Wise et al. 2018; Dumusque 2018; Chaplin et al. 2019; and references therein) and stellar activity has already led to several controversial planet detections in the past (e.g., αCenB b, Dumusque et al. 2012; Hatzes 2013; Rajpaul et al. 2016, GJ581 d and g, Vogt et al. 2010; Robertson et al. 2014, GJ667 c, and f Anglada-Escudé et al. 2013; Robertson & Mahadevan 2014).

This activity results from the contribution of various phenomena, which can be classified as a function of their correlation timescales. For main-sequence Solar-like stars, the three main noise sources originate from: (1) cyclic stellar oscillation eigenmodes (a few minutes), (2) stochastic surface convection motions (minutes–hours), and (3) (quasi) periodic stellar activity – spots, plages, flares – modulated with the stellar rotation or cycle (days–years). We note that even if active regions are more frequent at the maximum phases of the cycle, they can have lifetimes that are shorter than the rotation period (Saar & Donahue 1997).

In this paper, we aim to consider how the influence of the stochastic noise due to stellar convective motion could be counteracted. We ignore other sources of stellar RV variations (e.g., oscillations and active regions) and consider them as already corrected in our time series (e.g., through activity-sensitive lines, Baliunas et al. 1995; Wise et al. 2018 or dedicated filtering technique, Chaplin et al. 2019).

Convective noise can significantly alter the detection of exoplanets at the sub-m s⁻¹ level (Meunier et al. 2015, 2017b; Meunier & Lagrange 2019; Cegla 2019). The main technique proposed so far for mitigating its contribution in the RV series down to some tenth of cm s⁻¹ consists of averaging several (typically two or three) measurements of a target star during a night and separating them by at least two hours (Dumusque et al. 2011; Collier Cameron et al. 2019). However, the convection acts as a correlated noise over timescales longer than two-to-four hours (and even longer for supergranulation) and some correlations remain by using this observational procedure (Meunier et al. 2015). Moreover, this technique is performed at the cost of a small number of data points per night, leading to a critical lack of knowledge of the statistical properties for stellar activity as awhole. Other methods for dealing with convection noise consist of modeling the stellar activity as a correlated noise when fitting for RV planetary Keplerian signatures. Examples of common empirical models that we can find in the literature are red (i.e., power law) noise (Feroz & Hobson 2014), moving averaged noise (Tuomi et al. 2014), or Gaussian processes (Rajpaul et al. 2015). In practice, these empirical modelings should be used with caution as their results may lead to different conclusions, as shown in a recent RV challenge (Dumusque et al. 2017).

In this study, we question the reliability of traditional methods for determining the statistical significance of the detection in the presenceof correlated noise. This significance is based on the value of the false alarm probability (FAP) of statistical tests (see review in Khan et al. 2017). Traditionally, the FAP is derived under the assumption that the noise within the data (or the data residuals) is an uncorrelated white Gaussian noise (WGN). In this work, we propose a new method to access the significance of the detection of (quasi-)periodic signals in the presence of a correlated noise, providing that we can generate reliable (non-parametric) time series of this noise. We propose to use state-of-the-art 3D magneto-hydrodynamical (MHD) simulations of stellar surfaces to generate the noise series. We note that such simulations have already been investigated to determine the impact of convection on exoplanet detection (Cegla et al. 2013, 2018, 2019). Our analysis focuses on evaluating the reliability of MHD simulations in reproducing the time series of solar convective noise and on investigating the statistical benefit of using such simulated RV for deriving accurate FAP. In this paper, the benefit of using these simulations in the detection process is described for the case of regularly sampled time series. The case of an irregular sampling will be developed in a second paper², whereas the present analytical studies can nevertheless provide a useful proxy of the performance that can be expected in the case of an irregular sampling close to regular (e.g., one point per night at roughly the same hour).

This paper is divided into five sections. In Sect. 2, we evaluate the realism of the 3D MHD simulations. In Sect. 3, we use the standardized periodogram and present several detection tests to exploit this periodogram. In particular, we discuss some tests that are new to the exoplanet field: namely, the Higher-Criticism (Donoho & Jin 2004) and the Berk-Jones tests (Berk & Jones 1979). In Sect. 4, we perform a numerical study to investigate the benefit of our procedure and present our conclusions in Sects. 5 and 6.

2 Simulated solar granulation noise

As a preamble, we aim to test the realism of 3D MHD simulated RV time series of convective (granulation) noise and compare it to the RV time series obtained using the spectrophotometer Global Oscillation at Low Frequencies (GOLF).

2.1 Space measurements of RV solar convective noise

Measurements from spaceborne missions represent an excellent opportunity to validate the simulated velocities of solar convection. Indeed, they are not affected by the alternation of day and night or any ground-based follow-up problems (e.g., the influence of the Earth’s atmosphere) and provide regularly sampled time series at high cadence.

Since 1996, the GOLF spectrophotometer on board the Solar and Heliospheric Observatory (SoHO) spacecraft takes an almost continuous measurement of the solar disk-integrated position of the Sodium doublet. More particularly, it measures the position in the “blue” and “red” wings of the lines at ± 108 Å from the center of the lines, which are located at λ = 5895.924 Å ($D_1$) and 5889.950 Å ($D_2$). The solar light enters into a sodium vapor cell and a magnetic field splits the absorption lines (Zeeman effect). Then the Doppler shift (i.e., velocity) is evaluated as the flux ratio on these two points of the lines’ wings (see Unno et al. 1989, p. 328): $$v(t)\propto \frac{F_{\text{B}}(t)-F_{\text{R}}(t)}{F_{\text{B}}(t)+F_{\text{R}}(t)},$$(1)

where F_B and F_R are the fluxes inthe blue and red wings, respectively. For more technical details about this velocity extraction, we refer to Boumier & Dame (1993), Gabriel et al. (1995), Garcia et al. (2005) and Appourchaux et al. (2018).

After roughly one year of GOLF measurements, an instrumental failure happened and the velocity extraction was done using only one side of the sodium doublet: the blue wings (where the solar intensity comes from the bottom of the photosphere) from 1996 to 1998 and from 2002 until now and the red wings (where the solar intensity comes from the upper layers of the photosphere) between these dates (Garcia et al. 2005). Therefore, a careful calibration of the GOLF data was needed to obtain consistent velocities and several calibrations have been proposed. We chose to use the recent level-2 GOLF data³ calibrated as described in Appourchaux et al. (2018). In order to have the same sampling as in our MHD simulations, we sampled the GOLF time series every minute⁴. Moreover, we divided the GOLF time series into two-day sequences to study the RV correlations on the granulation timescales (from a few minutes up to several hours) and to validate them with Monte Carlo (MC) simulations based on a large number of solar subseries (see the statistical results presented in Sect. 3). From the entire sample of two-day sequences, we removed the ones containing observation gaps to have a perfectly regularly sampled time series (as our working hypothesis throughout this paper).

Finally, we applied a low-pass filter of 1620 μHz (i.e., 10.3 minutes) passband to filter out the oscillations modes and to restrict sensitivity to pick up only the convective noise. We computed the velocity root-mean-square (rms) of each 182 sequences available on the 1996 dataset (i.e., at solar cycle minimum, no calibration problem) and obtained an average value of 49 cm s⁻¹, which is in agreement with Pallé et al. (1999). An example of a two-day sequence and the corresponding periodogram will be shown in Sect. 2.2.4.

2.2 Synthetic time series of the RV solar convective noise

2.2.1 MHD simulations of the solar surface

We use the state-of-the-art radiative MHD code (stagger code, Nordlund & Galsgaard 1995) to simulate the surface convection and stratification of the Sun. In a 3D local-box model of the solar atmosphere (size: 8000 × 8000 kms and + 500 and − 3400 km above and below the surface at optical depth τ = 1), the code solves the full set of conservative MHD equations coupled to an accurate treatment of the radiative transfer. The horizontal sizes of the domain are defined to contain a sufficient number of granules at each time-step. The code is based on a sixth-order explicit finite difference scheme. The equations are solved on a staggered mesh where the thermodynamic variables are cell-centered, while the flux is shifted to the cell edge. The domain of simulation contains the entropy minimum located at the surface (photosphere) and is extended deep enough to have a flat entropy profile at the bottom (adiabatic regime). The code uses periodic boundary conditions horizontally and opened boundaries vertically. At the bottom of the simulation, the inflows haveconstant entropy and pressure. The outflows are not constrained and are free to pass through the boundary. We used a realistic equation-of-state that accounts for ionization, recombination, dissociation (Mihalas et al. 1988), and continuous line opacity (Gustafsson et al. 2008). Radiative transfer is crucial since it drives convection through entropy losses at the surface (Stein & Nordlund 1998) and is solved using the Feautrier’s scheme along with several inclined rays (one vertical, eight inclined) through each grid point. The wavelength dependence of the radiative transfer is taken into account using a binning scheme, in which the monochromatic lines are collected into 12 bins. The numerical resolution used for the present simulation is 120³. The choice of this modest resolution is a compromise between sufficient fine grid to catch enough of the inhomogeneities and sufficiently small to minimize the computing and storing costs of very long-run simulation. The solar parameters that define our 3D model are T_eff = 5775 ± 30 K, logg = 4.44 and a solar chemical composition (Asplund et al. 2009). The uncertainty in T_eff represents the fluctuations due to convection and p-modes. The average magnetic field in our simulation is ~100 G, as observed by (Hanle) spectropolarimetry (Trujillo Bueno et al. 2004).

In this work, we use an exceptionally long series of solar snapshots computed to study the properties of solar p-modes (Bigot et al., in prep.). It represents 53.14 days with a sampling of 60 s (i.e., 76 528 snapshots). To our knowledge, this is the longest series ever generated with such a 3D code. For the present study, we filtered out these modes since they have unrealistic large amplitudes (due to their small inertia) in such shallow boxes of granulation simulation. The synthetic sodium doublet lines are obtained for each snapshot by a monochromatic line transfer within [5884.000, 5901.945]  Å and at a resolution of20 000.

The synthetic line intensities I(t, x, y, λ, μ, ϕ) and the continuum C(t, x, y, λ, μ, ϕ) are computed for each x and y, the horizontal Cartesian coordinates of the simulation box and for several inclined rays defined by μ, the cosine of the six limb angles, and four azimuthal angles ϕ. The chosen discrete μ_i valuesare defined by the Gauss-Radau procedure. For six angles, we then have μ_i = {0.12, 0.39, 0.60, 0.80, 0.92, 1.00}. We averaged these intensities both horizontally and in azimuth to obtain our time-dependent center-to-limb intensity I(t, λ, μ) and continuum C(t, λ, μ), from which we will extract the radial velocities the following sub-sections.

Fig. 1 Comparison of velocities time series extracted from the Sodium doublet lines at different μ. For each time series, the oscillation modes have been filtered out.

2.2.2 RVdependence on the center-to-limb position

The radial velocities associated to each value of μ_i are obtained using Eq. (1). To compute Eq. (1), we generated the fluxes F(t, λ, μ_i) as the ratio of I(t, λ, μ_i) over C(t, λ, μ_i) for each μ_i. We then extracted the mean line profile F₀(λ, μ_i) by averaging the fluxes F(t, λ, μ_i) over t and used this reference profile to evaluate the fluxes ratio involved in Eq. (1). Finally, we translated these Doppler shifts into velocities using a proportionalfactor (κ) that results from a Taylor development around the considered wavelength λ₀ (see Unno et al. 1989, p. 328). This factor needs to be evaluated for each line of the Sodium doublet. It writes: $$\frac{1}{\kappa }=\frac{1}{c} \frac{\partial \text{ln}F(t,\lambda ,\mu )}{\partial \text{ln}\lambda }{|}_{\lambda ={\lambda }_0},$$(2)

with c the speed of light and λ₀ the wavelength corresponding to an intensity level of reference. We set this level of reference to $F_0({\lambda }_0^B,{\mu }_i)=F_0({\lambda }_0^R,{\mu }_i)=0.5$ with ${\lambda }_0^B$ and ${\lambda }_0^R$ the wavelengths in the blue and the red wings, respectively. The RV time series associated with three of the discrete {μ_i }-values are shown in Fig. 1.

The extracted radial velocities are strongly decreasing from the limb to the disk center, with typical rms velocities of 234.2 m s⁻¹ at the limb (μ = 0.12) and 23.3 m s⁻¹ at the disk center (μ = 1), as it is observed for the Sun (Löhner-Böttcher et al. 2018). This is explained by the fact that the observer does not see the same components of convective flows at the limb and the disk center. Indeed, the radial velocity is the projection of the total convective velocity, which includes both the vertical and horizontal velocities. At the disk center, radial and vertical velocities are the same, but at the limb, we only see the horizontal component. Since convection is strongly decelerating and horizontally diverging at the surface so that the gas overturns back to the interior in vertical downflows, the horizontal speeds are much larger than the vertical ones (Stein & Nordlund 1998; Nordlund et al. 2009). This explains the much larger values found at the limb than at the disk center. Moreover, the contribution of the small vertical velocity at the limb is strictly zero due to the projection effect. Despite the large rms velocities at the limb, we see in the following section that its contribution is limited due to surface projection effect when considering the disk-integrated velocities.

2.2.3 Disk-integrated RV

A single 3Dsimulation box represents a tiny fraction of the solar surface. Typically, we need $N_{\text{B}}=2\pi R_{\odot }^2/{\mathcal{l}}^2\approx 4.7\times {10}^4$ simulation boxes to cover the visible solar disk (i.e., half of the solar surface) with R_⊙ the solar radius and ℓ = 8 Mm the horizontal size of the simulation box. The difference in RV amplitude between those extracted from one simulation box (with μ-dependent rms velocity ≫ 1 m s⁻¹; see Fig. 1) and the solar disk-integrated observations (~49 cm s⁻¹) is due to the cancellation of positive (upflows) and negative (downflows) fluctuations when averaging over the entire disk. The reader might consult Ludwig (2006) for an in-depth discussion about this effect in the case of brightness fluctuations (see also Schrijver & Zwaan 2008). To generate a synthetic time series of the solar granulation as seen from disk-integrated observations, we follow a similar methodology to Ludwig (2006), Chiavassa et al. (2017) and Cegla et al. (2019) based on our single simulation box. The idea is to use the 76 528 available synthetic line profiles to patch a surface equivalent to the solar disk. Contrary to the previously mentioned studies, which used very short time series of a couple of hours, we have approximately 2 × N_B boxes to patch the solar disk at a given time, t, that allows us to cover the entire disk without duplication of the same snapshots. This allows us to avoid using the same snapshot in the patching procedure multiple times, which could lead to unavoidable correlations. In the present study, we randomly distributed the snapshots all over the surface, with the condition that two consecutive patches should correspond to times that are separated by at least 20 min to minimize possible correlations. For each patch k, we have an associated emergent intensity I(t, λ, μ_k). The values of μ_k are calculated depending to the position of the patch on the grid and the intensities I(t, λ, μ_k) interpolated from the six Gauss-Radau values using a second-order polynomial function. Then we let each of these patches evolves independently; that is, for each patch k, we performed a new interpolation from the six Gauss-Radau values to derive the new value for I(t + 1, λ, μ_k) corresponding to the considered μ_k. For each t and λ, we evaluated the disk integrated emergent flux as a function of wavelength, $$F(t,\lambda )={\sum }_{k=1}^{N_{\text{B}}}I(t,\lambda ,{\mu }_k) {\mu }_k,$$(3)

and we normalized the flux Eq. (3) by its corresponding value in the continuum $F_{\text{C}}(t,\lambda )={\sum }_{k=1}^{N_{\text{B}}}C(t,\lambda ,{\mu }_k) {\mu }_k$. As in Sect. 2.2.3, we then generated the mean line profile F₀ (λ), which is ourreference spectral line, to calculate the final Doppler shifts resulting from these disk-integrated synthetic Sodium line spectra. Finally, we extracted the Doppler velocity by measuring the flux ratio in the two points of each of the lines’ wings using Eq. (1) with the proportional factor given in Eq. (2).

The acoustic modes are naturally generated by the convective fluctuations inside a simulation box. However, in one shallow box, the modes have much lower inertia than the real observed p-modes. They have therefore much larger amplitudes. Hence, we eliminated their contribution to the RV time series by using a low-frequency filter of 1620 μHz passband (the same applied to GOLF time series; see Sect. 2.1). The rms of the final synthetic RV time series is 0.507 m s⁻¹, which is a value very close to the observed rms from space with GOLF (Pallé et al. 1999), and from the ground with HARPS-N (Collier Cameron et al. 2019). We note that our rms value does not take into account the possible contribution of the granulation noise to the high frequencies (ν > 1620 μHz) as we filtered them to remove the contribution of the stellar p-modes.

Other rms values due to granulation can be found in the literature. For example, Meunier et al. (2015) derived an rms that is twice higher (80 cm s⁻¹), Cegla et al. (2012) derived a similar value (40 cm s⁻¹) and Cegla et al. (2019) derived a smaller value (10 cm s⁻¹). The latter authors discuss the influence of the magnetic field that can reduce the velocity of the granulation flows in the 3D simulations.

2.2.4 Comparison between RV observations and simulations

To compare the synthetic velocity time series extracted from our 3D simulations with the GOLF observations, we added to both datasets a synthetic WGN to replace the high-frequency part of their power spectral density (PSD) that had been filtered out due to the presenceof the acoustic modes (see Sect. 2.2.3). The variance of this high-frequency noise was evaluated using the PSD at ν > 1620 μHz of the non-filtered GOLF data. This WGN does not affect the lower frequency part of the periodogram. We note that the influence of the four second exposure time of GOLF has been neglected in the computation of the PSD of the synthetic velocities.

The final comparison of the velocities is shown in the left panel of Fig. 2 for two selected two-day sequences. The corresponding (averaged) periodograms (see Eq. (8) in Sect. 3.2), resulting from the average of L = 26 regularly sampled two-day sequences, are shown in the right panel. The third periodogram represents the PSD of GOLF observations before the filtering of the acoustic modes. Toward the lower frequencies, we observe the frequency-dependent behavior of the solar granulation in all periodograms. When using high-resolution observations, the RV contribution of the stellar granulation acts as a frequency-dependent noise source that drastically differs from a WGN (characterized by a flat power over all frequencies). We observe a good match between the PSDs of the observed and simulated velocities until ν = 56 μHz (i.e., ~ 5 hr) corresponding to the correlation regime dominated by the granulation process (i.e., ν ∈ [50, 1000] μHz). At lower frequencies, the solar PSD becomes dominated by supergranulation and magnetic activity phenomena (spots and plages) and the comparison of these simulations of granulation with observations becomes obsolete (even if the granulation signal affects also the low frequencies of the PSD). We note that supergranules have longer lifetimes and should generate noise correlated over several days. They are not included in our present MHD simulations but they are also reproducible through 3D simulations, although their computing takes a longer time (e.g., Stein et al. 2009).

Fig. 2 Left: comparison of observed GOLF solar velocities (black) and synthetic velocities extracted from 3D simulationsof the granulation (red). The acoustics modes have been filtered out using a low-pass filter of 1620 μHz passband and a WGN has been added to both series. Right: associated averaged periodograms computed with L = 26 time series of 2 days duration. The grey PSD shows the averaged periodogram resulting from the unfiltered GOLF observations, where we can see the acoustics modes velocity signatures around 2000− 6000 μHz. The dotted line indicates the frequency regime where the high-frequency noise has been artificially added to both time series. The dashed line indicates the frequency limit (ν = 56 μHz) where the PSD is no longer dominated by the granulation noise.

3 Detection

This section presents the considered statistical model and detection tests. Several results detailed in Sulis et al. (2017b) are summarized below for the sake of completeness since Sect. 4 is aimed at validating the theoretical results from this study on real astrophysical data. The purpose of the approach is to detect (possibly quasi-) periodic components in a stationary colored noise with partially unknown statistics. By “partially” we mean that a training dataset of this colored noise is available (through MHD simulations). This noise dataset is independent of the observations (see Sect. 2). In the following, we assume the training data set contains all the noise sources that can affect the dataset under test. We note that currently, the MHD simulations per se cannot reflect the presence of active regions (spots, plages) due to the finite model precision. Hence, this study shows what can be done in the absence of such noise sources or in the situation where activity signatures can be identified by other means and added to the simulation.

3.1 Hypothesis testing problem

Let us consider a time series X(t_j) with N points, evenly sampled on times t_j = j ×dt for j = 1, …, N with dt the sampling time step. We consider a binary hypothesis problem of the form: $$H_0 :X(t_j)\text{ \&}=E(t_j),H_1 :X(t_j)\text{ \&}=R(t_j)+E(t_j),\text{.}$$(4)

where, under the null hypothesis, $H_0$, the data containonly the colored noise E(t_j) (of which a training set is available). The noise E is defined as a zero-mean second-order stationary Gaussian noise with unknown power spectral density S_E and absolutely integrable autocorrelation function r_E (see Sulis et al. 2017b). The alternative hypothesis, $H_1$, represents the case where an unknown RV planetary signal R(t_j) is melded with the colored noise. As illustrated, for example in Sulis et al. (2016), RV Keplerian signatures can be well approximated by a limited number of pure oscillations : $$R(t_j;{\theta }_{\text{R}})={\sum }_{q=1}^{N_{\text{s}}}{\alpha }_{\text{q}}\sin (2\pi f_{\text{q}}{t}_j+{\phi }_{\text{q}}),$$(5)

where the vector θ_R collects all the unknown amplitudes ${\alpha }_{\text{q}}\in {ℝ}^{+*}$, frequencies $f_{\text{q}}\in {ℝ}^{+*}$, and phases, φ_q ∈ [0, 2π[, of the N_s sinusoids. If a star reflects N_p planetary signatures, N_s is, in general, larger than N_p. The case N_s ≈ N_p corresponds to N_p planets with circular orbits and frequencies close to the Fourier grid. In all situations, N_s is much smaller than the number of Fourier frequencies.

3.2 Detection approach: a standardized periodogram

For simplicity, we consider for this section a unit time sampling d t = 1 and N even. When the observation sampling is regular, the search of periodic components can be done using the classical periodogram (Schuster 1898) defined as: $$P(\nu ):=\frac{1}{N} | {\sum }_{j=1}^{N}X(j) {\text{e}}^{-\text{i}2\pi \nu j}{|}^2.$$(6)

We note that to express P in units of density (m² s⁻² Hz⁻¹), expression (6) has to be divided by the passband 1∕d t (see Percival 1994, Eq. (11.6), Chap. 11). This leads us to consider in Eq. (6) a discrete Fourier frequencies defined as: $${\nu }_k:=\frac{k}{N}, \text{for} k=0,\dots ,N-1.$$

Owing to the hermitian symmetry of the Fourier transform and because we are not interested in the null frequency, below we consider P(ν) in Eq. (6) only as evaluated on a subset of $\frac{N}{2}-1$ independent Fourier frequencies corresponding to $k\in \Omega :=\{1,\dots ,\frac{N}{2}-1\}$. Asymptotically, P is an unbiased⁵ (but inconsistent⁶) estimate of the PSD (see Brillinger 1981, Thms. 5.2.1 and 5.2.4). The asymptotic distributions of P under both hypotheses are known ∀k  ∈ Ω: $$P({\nu }_k|H_0)~\frac{S_E({\nu }_k)}{2}{\chi }_2^2,(\text{see Brillinger 1981, Thm. 5.2.6)}),P({\nu }_k|H_1)~\frac{S_E({\nu }_k)}{2}{\chi }_{2,{\lambda }_k}^2,(\text{see Li 2014, Cor. 6.2}),$$(7)

with S_E the (unknown) noise PSD and λ_k = λ(ν_k;S_E, θ_R) a non-centrality parameter. For N_s sinusoidal components involved under $H_1$, the expression of this parameter can be found in Eq. (6) of Sulis et al. (2017b). We note that if the noise PSD S_E is unknown, the distribution of P given in Eq.(7) is also unknown.

Assuming now that L time series of the colored noise (denoted by {X_ℓ}, ℓ = 1, …, L below), can be generated under $H_0$ as a trainingdataset, we propose to use them as an estimate of the noise PSD to calibrate the periodogram of the data under test. Based on these L time series, we compute an averaged periodogram defined as: $${\overline{P}}_L({\nu }_k|H_0):=\frac{1}{L} {\sum }_{\mathcal{l}=1}^L \frac{1}{N} |{\sum }_{j=1}^N X_{\mathcal{l}}(j) {\text{e}}^{-\text{i}2\pi {\nu }_{k}j}{|}^2.$$(8)

We note that this periodogram has been initially introduced by Bartlett (1950), and used on subseries) to reduce the variance of the classical periodogram given in Eq. (6).

This averaged periodogram is an asymptotically consistent and unbiased estimator of the PSD. Following the same reasoning as for Eq. (7), the asymptotic distribution of ${\overline{P}}_L$ can be easily derived ∀k ∈ Ω as: $${\overline{P}}_L({\nu }_k|H_0)~S_E({\nu }_k) \frac{{\chi }_{2L}^2}{2L}.$$(9)

Using Eq. (8) to calibrate Eq. (6), we define the standardized periodogram as: $$\tilde{P}({\nu }_k):=\frac{P({\nu }_k)}{{\overline{P}}_L({\nu }_k)}.$$(10)

Thanks to the known distributions of the numerator and denominator of Eq. (10) and to their mutual independence, we can also derive the distribution of the standardized periodogram. Using Eqs. (7) and (9), we obtain a ratio of two independent χ² variables. This ratio leads, under $H_0$ and $H_1$, to a central and non-central F-distribution with respectively 2 and 2L degrees of freedom: $$\tilde{P}({\nu }_k|H_0)~\frac{{\chi }_2^2/2}{{\chi }_{2L}^2/2L}~F(2,2L),\tilde{P}({\nu }_k|H_1)~\frac{{\chi }_{2,{\lambda }_k}^2/2}{{\chi }_{2L}^2/2L}~F_{{\lambda }_k}(2,2L).$$(11)

We note that under ${\mathcal{H}}_0$, the distribution of the standardized periodogram $\tilde{P}$ is now asymptotically independent of the noise PSD S_E. This important property makes tests applied to $\tilde{P}$ act as Constant False Alarm Rate detectors (Scharf & Friedlander 1994): their false alarm rate is independent of the noise PSD. This is a very desirable feature in practice since it allows to control the false positive rate despite the unknown noise PSD. Under ${\mathcal{H}}_1$, the distribution depends on the noise PSD through the non-centrality parameter λ_k. The definition and the analysis of the theoretical performance of tests based on Eq. (10) are summarized in the following section.

3.3 Analysis of tests applied to the standardized periodogram

Before introducing the tests, it is convenient to consider vectors of random variables, noted in bold. For instance, the vector collecting the periodogram ordinates is written as: $$P:={[P({\nu }_1),\dots ,P({\nu }_N)]}^{\top }.$$

Notation x|y denotes a standardization of the entries of x by those of y. For instance, the vector of periodogram ordinates is standardized as in Eq. (10) and defined on the frequency set Ω. It is written as: $$\tilde{P}|{\overline{P}}_L:={\left[\frac{P({\nu }_1)}{{\overline{P}}_L({\nu }_1)},\dots ,\frac{P({\nu }_{\frac{N}{2}-1})}{{\overline{P}}_L({\nu }_{\frac{N}{2}-1})}\right]}^{\top }.$$

3.3.1 Test designed for a single periodicity

A common test consists of comparing the maximum periodogram value to a detection threshold $\gamma \in {ℝ}^{+}$ that determines the false alarm rate: $$T_M(\tilde{P}|{\overline{P}}_L):={\max }_k \tilde{P}({\nu }_k) \underset{{\mathcal{H}}_0}{\overset{{\mathcal{H}}_1}{\gtrless }}\gamma .$$(12)

This test is most efficient when a single periodicity on the Fourier grid is present under ${\mathcal{H}}_1$ (Donoho & Jin 2004). As the asymptotic distribution of $\tilde{P}$ is known at each frequency (see Eq. (11)), the false alarm and detection probabilities (noted P_FA and P_DET respectively), as well as their relationship (P_DET(P_FA)), can be derived analytically (Sulis et al. 2017b): $$P_{\text{FA}}(\gamma ):=\text{Pr}\left(T_{\text{M}}(\tilde{P}|{\overline{P}}_L)>\gamma |{\mathcal{H}}_0\right)=1-{\left(1-{\left(\frac{L}{\gamma +L}\right)}^L\right)}^{N_i},$$(13) $$P_{\text{DET}}(\gamma ):=\text{Pr}\left(T_{\text{M}}(\tilde{P}|{\overline{P}}_L)>\gamma |{\mathcal{H}}_1\right)\approx 1-{\prod }_{k\in \Omega }{\Phi }_{F_{{\lambda }_k}}(\gamma ,2,2L),$$(14) $$P_{\text{DET}}(P_{\text{FA}})\approx 1-{\prod }_{k\in \Omega }{\Phi }_{F_{{\lambda }_k}}(ϵ,2,2L),$$(15)

where $N_i:=\frac{N}{2}-1$ is the number of frequencies effectively considered in the test, $ϵ:=L[{(1-{(1-P_{\text{FA}})}^{\frac{1}{N_i}})}^{-\frac{1}{L}}-1]$ and ${\Phi }_{F_{{\lambda }_k}}$ is the cumulative distribution function (CDF) of a non-central F variable with non centrality parameter λ_k. The P_DET expressions given in Eqs. (14) and (15) are approximations due to the approximate independence of the periodogram ordinates under ${\mathcal{H}}_1$ (see Li 2014, Thm. 6.5). However, the analytic formulae above are quite accurate for values of N_i considered in practice as shown in Sulis et al. (2017b). These results allow saving a substantial amount of computation time for comparing the tests (in comparison with a MC simulation-based approach). They also allow gaining theoretical insight into the relative performances of the tests. Using the relation P_DET(P_FA), receiver operating characteristic (ROC) curves can be computed to compare the performances of the statistical tests. Furthermore, these analytical results can also be used to design detectability studies (see Sect. 4.2).

3.3.2 Tests designed for multiple periodicities

Testing for the largest peak in the periodogram may not be the best strategy for the case of multiple (quasi-) periodic signals. Chiu (1989) showed that in such cases tests exploiting order statistics of the periodogram may be more powerful than T_M (which looks at the maximum value only). In the case where the number of periodogram ordinates at Fourier frequencies affected by the planetary signature can be guessed or estimated a priori (let N_C denote this number), a generalization of test T_M replaces the maximum by the $N_{\text{C}}^{\text{th}}$ largest periodogram components. For such a test, analytic expressions for both the P_DET and P_FA can also be derived (see test T_C in Sulis et al. 2017b).

In practice, however, $N_{\text{C}}^{\text{th}}$ is often unknown and it is necessary to turn to tests that are adaptive with regard to the number of periodicities contained in the total Keplerian signature. Such tests are based on the P-values (noted v below) of the standardized periodogram. In the framework considered here, the P-values of an observed random variable (periodogram, or test statistic) is defined as the probability, under the null hypothesis, of obtaining a more extreme value than the observed one. Precisely, the P-values of $\tilde{P}({\nu }_k)$ are defined ∀k  ∈ Ω as: $$v_{\tilde{P},k}:=1-{\Phi }_F\left(\tilde{P}({\nu }_k),2,2L\right),$$

with Φ_F the CDF ofa central F variable. Examples of adaptive tests based on the P-values are the Higher-Criticism (Donoho & Jin 2004; Sulis et al. 2017b) and the Berk-Jones tests (Berk & Jones 1979; Aldor-Noiman et al. 2013; Mary & Ferrari 2014; Kaplan & Goldman 2014; Gontscharuk et al. 2015; Moscovich et al. 2016) respectively defined as: $$\text{HC}(\tilde{P}|{\overline{P}}_L):={\max }_{1\le k\le {\alpha }_{0}N}\frac{\sqrt{N}(k/N-v_{\tilde{P},(k)})}{\sqrt{v_{\tilde{P},(k)}(1-v_{\tilde{P},(k)})}} \underset{{\mathcal{H}}_0}{\overset{{\mathcal{H}}_1}{\gtrless }}\gamma ,$$(16)

and $$\text{BJ}(\tilde{P}|{\overline{P}}_L):={\max }_{1\le k\le {\alpha }_{0}N} I_{1-v_{\tilde{P},(k)}}(N-k+1,k)\underset{{\mathcal{H}}_0}{\overset{{\mathcal{H}}_1}{\gtrless }}\gamma ,$$(17)

where $v_{\tilde{P},(k)}$ denotes the order statistics of the P-values of the standardized periodogram (which have a beta distribution, David & Nagaraja 2003), α₀ is a constant $\in [\frac{1}{N},1],$ and I denotes the CDF of a beta variable.

Such tests as HC or BJ consist of setting a multiple testing problem, in which a set of test statistics (in our case, this refers to the periodogram at different frequencies) is taken and each of them are simultaneously considered in order to discriminate between the two hypotheses. In essence, these tests compare the maximal deviation of the empirical CDF of the ordered periodogram’s P-values to their true CDF under ${\mathcal{H}}_0$. The definition of the deviation depends on the test; both can be seen as variants of a generic divergence (Zhang et al. 2017).

In periodograms under ${\mathcal{H}}_1$, the planetary signature affects only a small fraction of the total number of ordinates; furthermore, this is by only a very smallamount, leading to a very difficult “needle in a haystack” detection problem. Donoho & Jin (2004) and Moscovich et al. (2016) demonstrate theoretically that HC and BJ present optimal guarantees in this regime. For finite values of N_i, the studies of Zhang et al. (2017) and Sulis et al. (2017b) show that BJ can be more powerful than other tests in case of weak and non extremely sparse signatures (e.g., multiplanetary systems of small planets with off-grid orbital frequencies and with high eccentricity orbits). We note that, in the case of irregular sampling – that will be the subject of a second paper, RV planet signatures can be much less sparse in the Fourier domain than for regular sampling owing to the sidelobes of the spectral window. Interestingly, efficient and accurate analytic calculations for the distribution of several adaptive tests, such as the HC and BJ under the null and the alternative hypotheses, have been recently included in Zhang et al. (2017). These features make adaptive tests particularly interesting for exoplanets detection, as illustrated in the numerical study below.

4 Numerical study

In this section, we first evaluate the validity of the statistical method presented in Sect. 3 using the solar observed and synthetic RV time series presented in Sect. 2. In a second step, we perform detectability studies for different planet signatures in the presence of solar convective noise by exploiting our analytical results. Finally, we compare the power of classical and adaptive detection tests for different Keplerian signatures.

4.1 Control of the false alarm: comparison of methods

The first part of this numerical study aims to compare the reliability of different false alarm probability estimates. We compare in particular bootstrap approaches to periodogram standardization (assuming a noise training data set is available). For the sake of concision, we focus on one test: the test of the maximum (see Eq. (12)). As for the considered dataset, we selected the regularly sampled two-day GOLF time series that are available for the first ten years of GOLF observations. In this sample, we removed sequences that are affected by strong outliers due to instrumental defects. This corresponds to a set of N_series = 1640 GOLF times series, with N = 2880 data points each. As described in Sect. 2, we filtered out the acoustic modes and added to each time series a WGN of standard deviation σ = 49 cm s⁻¹. This dataset represents our sample of solar observations under ${\mathcal{H}}_0$, as none of them contains any signs of the Solar System planets (the shortest period, of Mercury, is ≈ 88 days or 1.31 × 10⁻⁷ Hz) nor thestellar oscillations modes (affecting mostly the frequencies in the range 1–5 × 10⁻³ Hz) that have been filtered out. In the following, we will run tests on the frequency range that is dominated by the granulation noise: ν ∈ [50−8333] μHz.

As discussed in the introduction, a traditional approach in RV planet detection for evaluating FAP thresholds is based on bootstrap procedures. These methods assume that the observations (or their residuals in the case where some periodicities have been removed) contain only noise and that this noise is further uncorrelated with unknown variance. The noise statistics are estimated from the observations (see e.g., Jenkins et al. 2013, Hobson et al. 2018, Trifonov et al. 2018, Ment et al. 2018). The FAP is evaluated by estimating the distribution of the test statistic using fake data, typically obtained by shuffling the data.

Let us consider first the (scaled) max test, $$T_{\text{M}}(2\tilde{P}|{\sigma }^2):={\max }_k 2 \frac{P({\nu }_k)}{{\sigma }^2} \underset{{\mathcal{H}}_0}{\overset{{\mathcal{H}}_1}{\gtrless }}\gamma ,$$(18)

which, by definition, has FAP defined as $$P_{\text{FA}}(\gamma ;T_{\text{M}}(2 \tilde{P}|{\sigma }^2):=1-{\Phi }_{\text{M}}(\gamma ),$$(19)

where Φ_M is the CDF of T_M.

If the data contains a pure WGN of known variance σ², it can be shown (see e.g., Sulis 2017, Sect. 2.4.2) that $${\Phi }_{\text{M}}={(1-{\text{e}}^{-\gamma /2})}^{N_i},$$(20)

with N_i = N∕2 − 1 the number of considered (independent) periodogram components.

In the case where the variance σ² is unknown, the Max test, taking an estimate of the variance, ${\widehat{\sigma }}^2$, uses $T_M(2\tilde{P}|{\widehat{\sigma }}^2)$ as a test statistic. The bootstrap procedure consists in this case of estimating the variance and repeating the following steps: (i) shuffle the observed time series, (ii) compute the resulting periodogram on the new data set, and (iii) evaluate the test statistics Eq. (18) with ${\widehat{\sigma }}^2$ replacing σ². After generating a large number of realizations of test’s statistics, the FAP is derived as in Eq. (19), with the empirical distribution ${\widehat{\Phi }}_{\text{M}}$ replacing Φ_M.

This numerical procedure gives good results when the noise is white and Gaussian. This is illustrated in panel a of Fig. 3. This panel shows three FAP as a function of the detection threshold for the Max test. First, the blue line shows the FAP of test $T_M(2\tilde{P}|{\sigma }^2)$: this is the case for which σ² is known and the FAP is obtained using Eqs. (20) in (19). Second, since the bootstrap procedure describes above estimates σ² for the time series and follows steps (i)–(iii) above, the estimated function P_FA(γ) depends on the original data set used to generate the “fake” data set obtained by shuffling. One FAP estimate obtained for one particular data set, a WGN with standard deviation σ = 49 cm s⁻¹, is shown by the red curve, while the true FAP of $T_M(2\tilde{P}|{\widehat{\sigma }}^2)$ (evaluated on a WGN of variance σ² instead of ${\widehat{\sigma }}^2$), is shown by the dark green curve. Third, if we investigate the dependence of the FAP estimate with the original dataset (used to estimate ${\widehat{\sigma }}^2$), we obtain the orange curves of the panel a: here we show 100 curves corresponding to 100 different original data sets. We see from this panel that the bootstrap procedure is quite stable with respect to the considered dataset.

When the noise is colored, the data shuffling breaks the correlations present within the data and the situation changes. This is shown in panel b. The red curve shows the FAP of the test $T_M(2\tilde{P}|{\widehat{\sigma }}^2)$ estimated by bootstrap on one particular time series of estimated variance ${\widehat{\sigma }}^2$. The true FAP of this test as estimated from the remaining N_series − 1 is shown by the dark green curve. The orange (resp. light green) curves show the same as the red (resp. dark green) curve for all other times series. In contrast to the WGN case, the evaluation of the FAP is not robust nor reliable in this case. Hence, if noise correlations are ignored, a classical bootstrap procedure may severely underestimate the FAP and derive irrelevant thresholds.

This is further illustrated in panels d and e of Fig. 3. In all bottom panels, the estimated distributions of test T_M are shown in grey. The empirical thresholds corresponding to FA rates of 1 and 10% are represented by the vertical solid and dotted lines, respectively. Their numerical values can be read in Table 1. Panel d shows this distribution as obtained from one bootstrap procedure in the case of WGN. In this case, the thresholds estimated by bootstrap are close to the values they should have to ensure the target FAP. However, for the case of solar observations(panel b), these estimates are incorrect and lead to FAP that can be an order of magnitude larger than the target value. For instance, for one series (see red curve in panel b), the bootstrap procedure derives for a FAP of 1% a threshold value of γ = 23.65, whereas thisvalue is clearly underestimated: at this threshold, the true FAP, as estimated using all other time series, is in the range of [70.3%, 95.1%] (see green light curves).

To conclude this part of the analysis of results on GOLF data, we now turn to test T_M applied to the standardized periodogram (see Eq. (12)). In this case, the theoretical FAP is known (although the noise DSP is analytically unknown) and given by Eq. (13). To verify this expression, we standardize each of the periodograms of the GOLF sequences by the averaged periodogram computed using the L = 20 noise training datasets generated by the MHD simulations of the granulation (see Sect. 2). We then apply test $T_{\text{M}}(\tilde{P}|{\overline{P}}_L)$ and derive the associated P_FA as in Eq. (13). The results are shown in panels c and f of Fig. 3. This time, we observe in both cases a very good match between the theoretical FAP and the empirical values (see also last columns of Table 1).

We conclude the presentation of this first study with a short discussion. Of course, our point is not to show that the bootstrap is doomed to fail in case of colored noise; rather, it might possible to design bootstrap procedures that would take benefit from a training data set (as the approach of panel c does) or would use pre-whitening to obtain more robust FAP estimates than shown in panels b and e. Our point here is primarily to show that noise correlation caused by stellar convection severely impacts FAP estimates and that the proposed approach based on standardization achieves the desired robustness in estimating the FAP. These results validate the MHD simulation-based standardization approach for the control of the FAP and in particular the accuracy of the analytic calculations for test $T_{\text{M}}(\tilde{P}|{\overline{P}}_L)$ on real data. Since the principle of the approach based on accurate MHD simulations would be unchanged for a different spectral type, these results suggest that it can be used for detecting exoplanets orbiting any type of convective star.

Fig. 3 Illustration of the reliability of different FAP estimates of test TM depending on the noise characteristics under $H_0$ and on the technique involved. Top: FAP as a function of the detection threshold γ in the case where the data under $H_0$ is a WGN of standard deviation σ = 49 cm s−1 (first column) or the colored solar time series of same variance (last two columns). Panel a: the blue curve corresponds to the FAP of test TM with known variance σ2 (see Eq. (20)). The red curve corresponds to the FAP of test $T_{\text{M}}(2\tilde{P}|{\widehat{\sigma }}^2)$ estimated by bootstrap using one estimate of the variance from one time series, ${\widehat{\sigma }}^2$. The dark green curve represents the true FAP of $T_{\text{M}}(2\tilde{P}|{\widehat{\sigma }}^2)$. The curves in orange show 100 FAP estimates of the same test but obtained for 100 different estimates of ${\widehat{\sigma }}^2$. Panel b: the red curve shows the FAP of test $T_{\text{M}}(2\tilde{P}|{\widehat{\sigma }}^2)$, evaluated by bootstrap on one GOLF solar time series with estimated variance ${\widehat{\sigma }}^2$. The dark green curve shows the true FAP of this test, as estimated using the Nseries− 1 other GOLF time series. The orange (resp. light green) curves are the same as the red (resp. the dark green) curves, but using each time a different GOLF time series as input. Panel c: the green curve represents the analytic FAP of TM based on the simulation-standardized periodogram with L = 20 MHD time series (see Eq.(13)) and the red curve represents the true FAP estimated using Nseries = 1640 GOLF time series. Bottom: empirical distribution of test statistics TM as estimated by bootstrap (panels d and e) and by MC simulations of the GOLF series standardized by the MHD simulations (panel f). In all six panels, the thresholds inferred for FAPs of 1 and 10% by each technique are indicated by the dashed and dotted lines, respectively. The color used for the thresholds in each bottom panel corresponds to the color used for each method in the corresponding upper panel. Numerical values are indicated in Table. 1.

4.2 Detectability study

As we have seen in Sect. 3, exploiting MHD simulations of the granulation noise opens up the possibility for analytically controlling the false alarm rate and extending the power of the tests for any values of the observation parameters. Comparing the impact of these parameters on the probability of detection for a fixed FAP is very useful in designing observational strategies, for instance.

Let us consider again the test T_M given in Eq. (12), for which the detection probability can be computed using expression (15) for a given P_FA. Given a specific planetary signature, we want to evaluate the observation duration (T_obs) that is required to allow for the detection of this planet with a large probability (say, P_DET = 80% at P_FA = 1%). We simulated for this study different planetary signatures under $H_1$ with circular orbits and orbital frequencies on the Fourier grid (we slightly adjusted the time sampling step d t as T_obs increases to guarantee that the period is exactly on the grid). For such signatures, only one periodogram ordinate is affected under ${\mathcal{H}}_1$, while T_M is optimal (Donoho & Jin 2004). For periodogram standardization, again we used the simulated velocities discussed in Sect. 2. The considered convection noise corresponds to a Sun-like star.

Some results are shown in Fig. 4. Each panel of the figure investigates the influence of a different parameter (see legend). The black curves correspond to a configuration where a 0.5 Earth-mass planet orbits circularly its host star with a period of 17.5 hr, the regular time sampling step dt is 2 hr and L = 20 MHD simulations time series are available for periodogram standardization. This setting corresponds to a RV signature of semi-amplitude K = 0.35 m s⁻¹ and an orbital frequency of f_p = 1.58 × 10⁻⁵ Hz. In each panel, the pale blue dot indicates a detection probability of 80% in this configuration for this planet. The dashed lines in the first three panels represent the detectability in the case where the noise is white (instead of colored) but with the same standard deviation as the colored convection noise (σ = 49 cm s⁻¹). The analytical expression for this probability is (see Eq. (2.53) in Sulis 2017): $$P_{\text{DET}}(\gamma ;T_M(2\tilde{P}|{\sigma }^2)):=1-{\prod }_{k\in \Omega }{\Phi }_{{\chi }_2^2,{\lambda }_k}(\gamma ),$$(21)

where ${\Phi }_{{\chi }_2^2,{\lambda }_k}$ is the CDF of a non central ${\chi }_2^2$ distribution with two degrees of freedom and the non-centrality parameter λ_k (see Eq. (6) of Sulis et al. 2017b).

The panel a of Fig. 4 shows, for instance, that in the considered configuration, an observation run totaling T_obs ≈ 12.4 days (corresponding to N ≈ 150 sample points with d t = 2 hr) would allow the detection of a 0.5 Earth-mass planet with a probability of 80%, while ensuring a false alarm rate of 1%. We note that, in contrast, we would only need T_obs = 3.0 days to reach the same trade-off P_DET vs. P_FA if the noise was uncorrelated (see the dashed circle). This factor ≈4 in duration is the price that has to be paid in order to fight against correlation caused by convection noise. If the planet mass is lower (panel a), the required observational time T_obs can increase extensively. For example, for a planet with a mass similar to Mars (≈ 0.1 M_⊕ leading to an RV semi-amplitude K = 0.07 m s⁻¹), we would need at least 457 days of observations to achieve the same performances. Similarly, if the planet’s period increases, the needed observational time increases (because the amplitude of the Keplerian signature decreases, which is not shown).

The test’s performance depends also on the sampling of the orbital frequency (panel b). In our example, if the orbital frequency is not on the Fourier frequency grid, the detection performance of this test decreases, with a loss in P_DET that can reach a factor of 2.

Increasing the sampling time step (panel c) or decreasing the number of used training data set (panel d), increases also T_obs. We also note that there is a very small improvement of the test performance brought by increasing L as soon as L is sufficiently large (for L = 50 and L = 1000, when the required observation durations are T_obs = 10.5 and 9.5 days, respectively). This fact is particularly interesting since the MHD simulations are computationally heavy and L = 1000 may remain outside of the reach of the coming decades.

These plots are examples of false alarm versus power trade-offs that can be achieved by exploiting reliable time series of the convective colored noise. We note that the values indicated in this study are drastically different from those reported in Sulis et al. (2017b). For instance, we reported T_obs = 250 days for an 1.1 M_⊕ planet orbiting its star in 3.2 days with dt = 4 h and L = 100, while with these same parameters, we find now T_obs ≈ 17 days. The reason is that the PSD considered to represent the solar granulation noise source is different from that given in these first works: the considered PSD is now more realistic and deeply checked against Solar observations (see Sect. 2).

We now give an example of an application of adaptive tests, which are less well known in the exoplanet community than test T_M, although they can sometimes present advantages over the latter.

Fig. 4 Detection probability as a function of the observation time for a single planet in circular orbit with period 17.5 hr around a solar-type star, for test $T_{\text{M}}(\tilde{P}|{\overline{P}}_L)$, at PFA= 1%. The orbital inclination is set to 90 degrees. The different panels show the influence of the planet’s mass (a), the orbital period (with the corresponding orbital frequency on (black) or off (gray) -Fourier grid) (b), the time sampling step (c) and the size of the training data set (i.e., the number of available noise times series) (d). In each panel, the black curve indicates the PDET obtained for a configuration in which the planet has a mass of 0.5 M⊕, a circular orbit and L = 20 HD time series are available for periodogram standardization. The point where PDET reaches 80% for this configuration is indicated by the blue disks. In each panel, the legends indicate the parameters under study. The dashed lines represent the planet’s detectability in the case where the noise is white (instead of colored) but with the same standard deviation as the colored convection noise (σ = 49 cm s−1).

4.3 Detectability of general Keplerian signatures

In this section, we compare the performances of the different tests presented in Sect. 3, i.e., T_M Eq. (12), HC Eq. (16) and BJ Eq. (17) for different types of Keplerian signatures. The combination of Keplerian parameters influences the shape of the RV signature which, in turn, influences the sparsity of the signature in the Fourier domain, that is, the number of periodogram components affected by the presence of a planetary signature (e.g., see Sulis et al. 2016 for a detailed study of the influence of Keplerian parameters on sparsity). Here we define the sparsity coefficient $S_{\beta }$ as the proportion of non-zero coefficients and, as in Donoho & Jin (2004), we parameterize $S_{\beta }$ as: $$S_{\beta }:=\frac{N_s}{N}:=N^{-\beta },$$

with β ∈ [0, 1] a sparsity parameter. The value β = 1 corresponds to an extremely sparse signal (i.e., a single periodogram frequency is affected by the periodic signal), and β → 0 to a non sparse signature. RV signatures correspond in general to sparse signatures (β is typically in the range $[\frac{1}{2},1]$). The less sparse signatures are obtained for multiple systems, planets having highly eccentric orbits and planets with off-Fourier grid orbital frequencies.

4.3.1 Adaptive tests

We compare the detection probability of tests T_M, HC and BJ for two types of planet signatures. For the first case, we consider the signal of a 0.4 M_⊕ planet in a circular orbit with frequency on the Fourier grid. In the second case, we consider the same planet but with a highly eccentric orbit (e = 0.9) and a slightly off-grid orbital frequency. The RV signatures of these two planets and the corresponding periodograms are shown in the top panels of Fig. 5. We see a difference in the number of significant peaks in their corresponding periodograms.

We performed MC simulations to evaluate the ROC curves of tests T_M, HC and BJ. As the number of available MHD simulations is limited, we generated the noise under $H_0$ using a 20 order autoregressive (AR) process, with parameters fitted to the MHD simulated time series. We generated 10⁴ realizations of this AR process under $H_0$ and 10⁴ other realizations under $H_1$ with the two RV planetary signatures. For each realization, we computed the standardized periodogram in Eq. (10) using L = 20 series for the denominator in Eq. (8), and applied the different tests. Results are shown in the bottom panel of Fig. 5. Comparing the performances of test T_M with the adaptive tests, we observe for the “sparse signal” the best results for T_M over the other two tests, with HC close to T_M at low FAP. However, in the case of a high-eccentricity planetary orbit, tests BJ and HC show better performances than T_M at all FAP. We see that these adaptive tests present another important side advantage over T_M : their test statistics can be used to reliably estimate the frequency content for complex planetary signatures.

Fig. 5 Top: synthetic RV time series (left) and corresponding periodograms (right) for a single planet in circular (black) or eccentric (red) orbit around a Solar-type star. For both signals, the planet mass is set to 0.4 M⊕, the period to17.5 hr, the orbital inclination to 90 degrees, the argument at periastron to π∕2 radian, the time sampling step to 2 hr and the number of data points is N = 300. For the eccentric planet, the planet orbital frequency is slightly off the Fourier-frequency grid. Bottom: ROC curves of tests TM (solid), HC (dashed) and BJ (dotted) applied to $\tilde{P}|{\overline{P}}_L$ with L = 20 for the considered circular (left) and eccentric (right) orbital signals.

4.3.2 Ability of adaptive tests in recovering the signal’s frequency support

Here we evaluate the ability of each test in detecting the correct number of periodic components at their true location (i.e., the true signal frequencies in the periodogram). We note that this problem is different from that of discriminating between the “noise only” vs. “planet plus noise” hypotheses. To be able to quantify easily the number of periodic components to be detected under $H_1$, we consider periodic signals in the form of a sum of N_s pure sinusoidal signals with frequencies on the Fourier frequency grid. Moreover, to generate a large amount of MC simulations, the colored noise is generated as a low order AR process. The detection thresholds for a target FAP of 5% were derived for all considered tests by MC simulations on 10⁴ noise sequences under $H_0$. For the training dataset used for periodogram standardization (see Eq. (10)), we generated L = 20 synthetic separate noise time series for each of these 10⁴ sequences.

In each case, we computed tests T_M Eq. (12), HC Eq. (16), and BJ Eq. (17) on $\tilde{P}|{\overline{P}}_L$. By definition, test T_M only focuses on the largest component of the (standardized) periodogram. In contrast, tests HC and BJ focus on one particular ordered P-value: the one for which the corresponding test statistic is maximum. This ordered P-value, say $v_{\tilde{P},i^{\star }}$ corresponds to the i^⋆ largest periodogram components. Consequently, when a detection is made, the i^⋆ largest periodogram ordinates can be used to estimate the signal’s frequency support. In the following, we estimate by MC simulations the probability that i^⋆ = N_s and that the identified frequencies correspond to the true signal frequencies, at a fixed FAP of 5%.

Under $H_1$, we varied the number N_s of sinusoidal signals (by varying parameter β) added to each of the colored noise time series: N_s ∈ [1, 100], corresponding to β ∈ [0.33, 1]. The N_s sinusoids’ amplitudes denoted by A in Fig. 6 were taken as equal, with A varied in the range [0.1, 4.5] m s⁻¹, while their N_s frequency locations f_q were picked randomly in the Fourier grid.

The results are shown in Fig. 6 as a function of the signal parameters (amplitude and sparsity) for tests T_M, HC, and BJ (from left to right, respectively). For T_M, we only show the case of β = 1 as this test statistic focuses on the largest periodogram component. In the left panel, test T_M appears to be the best to correctly locate the signal frequency when only one frequency is present (the adaptive tests are shown for comparison). For N_s = 1 (β = 1), the probability of correctly locating the sinusoid frequency (P_R) grows faster to 1 for this test than for the other tests (with close to T_M). In contrast to T_M, tests HC and BJ allow both for the detection and characterization of the planetary frequency support (middle and right panel of Fig. 6). We also note the particularly good performances of BJ vs. HC for recovering the frequency support over a large sparsity regime. These tests, which benefit from theoretical optimal results (see Donoho & Jin 2004 and Moscovich et al. 2016), can be exploitable in the context of exoplanet detection by RV thanks to the considered periodogram standardization. Their good performances in theory and practice make them particularly interesting for the detection and characterization of extrasolar planetary systems.

Fig. 6 Probability of detecting the correct support (i.e., to find the true number of sinusoids with the correct frequencies) for tests TM (left), HC (middle) and BJ (right) applied to $\tilde{P}|{\overline{P}}_L$ with L = 20. For each given value of β, the Ns amplitudesare equal. Left panel: represents this probability as a function of the sinusoid’s amplitude (Ns = β = 1) and the performances of tests HC and BJ have been added for comparison. The two other panels represent the adaptive tests’ performances as a function of the sinusoids’ amplitude A and the sparsity parameter β. In these twopanels, the probability of correct recovery is indicated in color.

5 Discussion

5.1 Scope of the proposed method and a zoom on USP planets

Convective noise affects all components of the periodogram but its effects impact mostly the frequency range corresponding to periods between some minutes and several hours for solar-like stars. In Sect. 2, we obtained a good match between the observed and simulated PSD of solar RV in this frequency range. In practice, the convective noise cannot be “corrected” as magnetic activity may be (e.g., with chromospheric indicators, Baliunas et al. 1995; Wise et al. 2018) and constitutes a noise barrier, for which the statistical properties need to be known to reliably claim any planet detection at the cm s⁻¹ level. This was the purpose of our study and its presentation of the formalism of the approach. In a subsequent work, we will apply this formalism to other Solar-like stars having different convective properties.

The proposed method of basing such studies on a standardized periodogram could be directly applied to improve the determination of the FAP in the case of ultra-short period (USP) planets (defined with periods < 1 day) under $H_1$. USP planets are known to be tidally locked to their host star (leading to circular orbit) and of small size or mass (< 10 M_⊕). They exist, in general, in multi-planetary systems. According to Winn et al. (2018), this category of planets is as frequent as hot-Jupiters (defined with periods ranging up to 10 days), with one over 200 Sun-like stars hosting such planets. To illustrate the performance of the proposed technique for detecting USP planets, we performed a similar detectability study as described in Sect. 4.2 for some known USP planets (we assume their hosting star is similar to the Sun). The results for a sampling rate of d t = 12 hr and for L =2, available synthetic noise light curves are reported in Table. 2 and displayed in Fig. 7. The figure represents the observational time required to reach a P_DET of 80% for a P_FA equal to 1%, using test T_M given in Eq. (12), as a function of the planet’s orbital period. The colored curves show the observational time for virtual planets of different periods and masses. The black crosses correspond to real planets. Logically, we observe the increase of the observational time with the decrease of the planet mass and the increase of the planet period. The table indicates that most USP planets are detectable at the levels specified above with our technique within a couple of weeks (< 21 days). For small mass USP planets, with M_p < 2 M_⊕, it would take a couple of months to achieve the same performances.

Finally, we expect the method presented here to be easily extended to larger period ranges by computing MHD supergranulation instead of granulation, that is, by extending the simulation domain (making it larger and deeper) with exactly the same simulation setup (Rincon & Rieutord 2018). Since this would be more demanding in terms of CPU and storage, while retaining, in principle, what is shown here, we restrict the scope of this study to granulation scale simulations.

Fig. 7 Observation time as a function of the planet orbital period that is needed to achieve PDET = 80% and PFA = 1% with test TM (colored curves). The fixed parameters are dt = 12 hr and L = 2. The observation time associated with some known USP exoplanets (see Table. 2, assuming a Solar-like star) are represented by black crosses.

5.2 Benefits and limitations of the method

MHD simulations are non parametric, meaning that they do not rely on any adjustable parameter to fit the observed data (see Sect. 2.1). In practice, the frequency dependence of the granulation is often estimated using parametric laws, such as Harvey-like profiles (Harvey 1985). However, as demonstrated in Sects. VI and VII.D of Sulis et al. (2017b), the estimation of the parameters of these models leads to the injection of an estimation noise in the detection process, the statistics of which are difficult to capture. Besides, the noise parameters derived in this way may be contaminated by the signal to be detected and the choice of the noise parametric model can be subjective. As we show in this study, using an MHD simulation-based approach allows us to accurately control the estimation noise through the number L, whose impact on the tests performances can be exhibited analytically (see Eqs. (14)–(13)).

This study needs to be extended to other convective stars. The impact of granulation changes throughout the HR diagram: the larger is the pressure scale height at the surface (that is for larger effective temperatures or lower gravities), the larger are the fluctuations induced by the convective motions. The current limitation in the present method is the computational cost of the MHD simulations to generate a substantially long time series of velocities. However, in the coming years, the increased speed of CPU resources will allow for such computations to be carried out in a more systematic way. In a subsequent work, we will explore these effectsfor selected targets in the HR diagram. Based on the realism of the 3D MHD simulations, our results suggest thatthe proposed method can be a powerful and reliable way of detecting RV exoplanet signatures at the cm s⁻¹ level in the presence of convective noise.

A second limitation is the regular sampling involved in this study. In practice, the RV data are irregularly sampled and the FAP of any test based on any periodogram (Schuster 1898; Scargle 1982; Zechmeister & Kürster 2009) cannot be controlled by analytical expressions in the case of correlated noise because the periodogram components are interdependent.

However, as mentioned earlier in this paper, we underline that if the irregularity of the considered sampling remains weak, the analytical studies presented here may provide a useful proxy of the tests’ performance in practical situations. For example, this can be used to design detectability studies. For strongly irregular samplings, the techniques based on the MHD standardized periodogram presented here need to be adapted by dedicated bootstrap procedures (see, i.e., Sulis et al. 2017a). The application of this procedure to the real data deserves a full study that will be the purpose of a second paper.

6 Conclusions

In cases where the effective temperature, surface gravity, and metallicity of the star are precisely known (thanks to asteroseismology, interferometry, or spectroscopy), 3D MHD simulations are capable of generating realistic RV time series of the stellar granulation. This has been demonstrated for the Sun as part of studies involving the comparison of velocities extracted from 3D spectra of the sodium doublet and GOLF/SoHO observations.

Following the theoretical analysis described in Sulis et al. (2017b), we used these synthetic time series of the granulation colored noise to design standardized periodograms. These new standardized periodograms allow for the application of tests that are both powerful and for which we can derive accurate FAP. We present extensive numerical results based on real and synthetic data, including studies on the robustness, the detectability, and the frequency support recovery. In particular, we introduced adaptive tests, which are new in the field of RV planet detection. Even if the objective of this study is to detect planets down to the cm s⁻¹ level, the proposed procedure is a general approach that can be applied to many periodicity detection problems in astrophysics (and beyond).

Acknowledgements

This work was supported by the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES. S. Sulis acknowledges support from the Austrian Research Promotion Agency (FFG) under project 859724 “GRAPPA”, as well as Thales Alenia Space and PACA region. D.M. acknowledges support from the GDR ISIS through the Projet exploratoire TASTY. Computations have been done on the “Mesocentre SIGAMM” machine, hosted by Observatoire de la Côte d’Azur. The GOLF instrument onboard SOHO is a cooperative effort of scientists, engineers, and technicians, to whom we are indebted. SOHO is a project of international collaboration between ESA and NASA.

References

All Tables

All Figures

	Fig. 1 Comparison of velocities time series extracted from the Sodium doublet lines at different μ. For each time series, the oscillation modes have been filtered out.
In the text

Fig. 2 Left: comparison of observed GOLF solar velocities (black) and synthetic velocities extracted from 3D simulationsof the granulation (red). The acoustics modes have been filtered out using a low-pass filter of 1620 μHz passband and a WGN has been added to both series. Right: associated averaged periodograms computed with L = 26 time series of 2 days duration. The grey PSD shows the averaged periodogram resulting from the unfiltered GOLF observations, where we can see the acoustics modes velocity signatures around 2000− 6000 μHz. The dotted line indicates the frequency regime where the high-frequency noise has been artificially added to both time series. The dashed line indicates the frequency limit (ν = 56 μHz) where the PSD is no longer dominated by the granulation noise.

In the text

Fig. 3 Illustration of the reliability of different FAP estimates of test TM depending on the noise characteristics under $H_0$ and on the technique involved. Top: FAP as a function of the detection threshold γ in the case where the data under $H_0$ is a WGN of standard deviation σ = 49 cm s−1 (first column) or the colored solar time series of same variance (last two columns). Panel a: the blue curve corresponds to the FAP of test TM with known variance σ2 (see Eq. (20)). The red curve corresponds to the FAP of test $T_{\text{M}}(2\tilde{P}|{\widehat{\sigma }}^2)$ estimated by bootstrap using one estimate of the variance from one time series, ${\widehat{\sigma }}^2$. The dark green curve represents the true FAP of $T_{\text{M}}(2\tilde{P}|{\widehat{\sigma }}^2)$. The curves in orange show 100 FAP estimates of the same test but obtained for 100 different estimates of ${\widehat{\sigma }}^2$. Panel b: the red curve shows the FAP of test $T_{\text{M}}(2\tilde{P}|{\widehat{\sigma }}^2)$, evaluated by bootstrap on one GOLF solar time series with estimated variance ${\widehat{\sigma }}^2$. The dark green curve shows the true FAP of this test, as estimated using the Nseries− 1 other GOLF time series. The orange (resp. light green) curves are the same as the red (resp. the dark green) curves, but using each time a different GOLF time series as input. Panel c: the green curve represents the analytic FAP of TM based on the simulation-standardized periodogram with L = 20 MHD time series (see Eq.(13)) and the red curve represents the true FAP estimated using Nseries = 1640 GOLF time series. Bottom: empirical distribution of test statistics TM as estimated by bootstrap (panels d and e) and by MC simulations of the GOLF series standardized by the MHD simulations (panel f). In all six panels, the thresholds inferred for FAPs of 1 and 10% by each technique are indicated by the dashed and dotted lines, respectively. The color used for the thresholds in each bottom panel corresponds to the color used for each method in the corresponding upper panel. Numerical values are indicated in Table. 1.

In the text

Fig. 4 Detection probability as a function of the observation time for a single planet in circular orbit with period 17.5 hr around a solar-type star, for test $T_{\text{M}}(\tilde{P}|{\overline{P}}_L)$, at PFA= 1%. The orbital inclination is set to 90 degrees. The different panels show the influence of the planet’s mass (a), the orbital period (with the corresponding orbital frequency on (black) or off (gray) -Fourier grid) (b), the time sampling step (c) and the size of the training data set (i.e., the number of available noise times series) (d). In each panel, the black curve indicates the PDET obtained for a configuration in which the planet has a mass of 0.5 M⊕, a circular orbit and L = 20 HD time series are available for periodogram standardization. The point where PDET reaches 80% for this configuration is indicated by the blue disks. In each panel, the legends indicate the parameters under study. The dashed lines represent the planet’s detectability in the case where the noise is white (instead of colored) but with the same standard deviation as the colored convection noise (σ = 49 cm s−1).

In the text

Fig. 5 Top: synthetic RV time series (left) and corresponding periodograms (right) for a single planet in circular (black) or eccentric (red) orbit around a Solar-type star. For both signals, the planet mass is set to 0.4 M⊕, the period to17.5 hr, the orbital inclination to 90 degrees, the argument at periastron to π∕2 radian, the time sampling step to 2 hr and the number of data points is N = 300. For the eccentric planet, the planet orbital frequency is slightly off the Fourier-frequency grid. Bottom: ROC curves of tests TM (solid), HC (dashed) and BJ (dotted) applied to $\tilde{P}|{\overline{P}}_L$ with L = 20 for the considered circular (left) and eccentric (right) orbital signals.

In the text

Fig. 6 Probability of detecting the correct support (i.e., to find the true number of sinusoids with the correct frequencies) for tests TM (left), HC (middle) and BJ (right) applied to $\tilde{P}|{\overline{P}}_L$ with L = 20. For each given value of β, the Ns amplitudesare equal. Left panel: represents this probability as a function of the sinusoid’s amplitude (Ns = β = 1) and the performances of tests HC and BJ have been added for comparison. The two other panels represent the adaptive tests’ performances as a function of the sinusoids’ amplitude A and the sparsity parameter β. In these twopanels, the probability of correct recovery is indicated in color.

In the text

	Fig. 7 Observation time as a function of the planet orbital period that is needed to achieve PDET = 80% and PFA = 1% with test TM (colored curves). The fixed parameters are dt = 12 hr and L = 2. The observation time associated with some known USP exoplanets (see Table. 2, assuming a Solar-like star) are represented by black crosses.
In the text