© ESO 2020

1. Introduction

Detecting periodic signals in unevenly spaced time series is a common problem in astronomy, which is encountered, for instance, when searching for binaries or exoplanet companions in radial velocity, astrometric, or photometric time series. The Lomb-Scargle (LS) periodogram (Lomb 1976; Scargle 1982) is a classical and efficient tool to search for sinusoidal signals. The principle of the LS periodogram is to scan a wide range of frequencies and to compare a linear sinusoidal model at a given frequency with a constant model, called the base model. A widespread variant of the LS periodogram, called the generalized LS periodogram (GLS), was proposed by Zechmeister & Kürster (2009) where the constant is adjusted for at each frequency (see also Ferraz-Mello 1981). Once the periodogram is computed, the false alarm probability (FAP) criterion is often used to determine whether or not it supports the detection of a periodic signal. The FAP is often estimated using numerical methods such as bootstrapping or Monte Carlo. These methods can be computationally intensive, especially at low FAP levels. Indeed, estimating a FAP level of P numerically requires at least the computation of 10/P periodograms.

Several analytical formula have been proposed for the FAP (Scargle 1982; Horne & Baliunas 1986), but have been subsequently contested (Koen 1990). The periodogram framework was generalized in a series of works to handle more complex models (Baluev 2008, 2009, 2013a,b, 2015), where rigorous and sharp analytical approximations of the FAP were provided based on the so-called Rice formula and previous works by Davies (1977, 1987, 2002). These works allow a fast and rigorous estimation of the periodogram FAP in the presence of white noise. However, the white noise assumption is often incorrect for astronomical time series. Indeed, several sources of correlated noise, such as the astronomical source itself, Earth’s atmosphere, or instrumental systematics, could contaminate the measurements.

For instance, stellar variability has a huge impact on the detection of low-mass exoplanets using high-precision radial velocity time series. Indeed, at a precision better than 1 m s⁻¹, stellar variability affects the measurements on timescales ranging from a few minutes (for stellar oscillation p-modes), to hours and days (for stellar granulation and super-granulation), and even up to the star rotation period (for the effect of spots and plages). In this context, low-mass exoplanet detection becomes a challenge (e.g., Queloz et al. 2001) and proper tools have to be developed to treat correlated noise properly and to compute reliable periodograms and FAPs. Sulis et al. (2016) provide an analytical FAP estimate for periodograms normalized by the power density spectrum of the noise, in the limit of low aliasing and in the case of evenly sampled data. This is, however, not the case of most radial velocity datasets. Baluev (2013c) provides a “suggestive generalization” to the correlated noise case of the FAP formula obtained in Baluev (2008), but advocates against the use of this formula since it has not been proved rigorously.

In this article, we extend the work of Baluev (2008) to account for correlated noise. In Sect. 2 we define a general periodogram for an arbitrary covariance of the noise and provide an analytical approximation of the corresponding FAP, which we validate against Monte Carlo simulations. In Sect. 3 we provide a method to explore the sensitivity of the periodogram to the noise model. We discuss our results in Sect. 4.

2. Significance of periodogram peaks with correlated noise

In this section, we present a method to assess the significance of periodogram peaks (FAP) in the correlated noise case. We first give a general definition of the periodogram in Sect. 2.1. We then provide formulas to compute the corresponding FAP in Sect. 2.2. Finally, we compare this analytical FAP with the results of Monte Carlo simulations in Sect. 2.3.

2.1. General linear periodogram

We extend the general definition of least squares periodograms by Baluev (2008) to the correlated noise case. Following Baluev (2008), we compare the χ² of the residuals of a linear base model ℋ of p parameters with enlarged linear models 𝒦 of p + d parameters, parameterized by the frequency ν. The base model ℋ is written as

$$\begin{array}{c}\hfill \mathcal{H}:m_{\mathcal{H}}({\theta }_{\mathcal{H}})={\phi }_{\mathcal{H}}{\theta }_{\mathcal{H}},\end{array}$$(1)

where θ_ℋ is the vector of size p of the model parameters, φ_ℋ is a n × p matrix, and n the number of points in the time series. The columns of φ_H are thus explanatory time series that are scaled by the linear parameters θ_ℋ. For instance, in the case of a radial velocity time series with two different instruments and a linear drift, the linear base model could be chosen as

$$\begin{array}{c}\hfill m_{\mathcal{H}}={\gamma }_1{\delta }_1(t)+{\gamma }_2{\delta }_2(t)+\alpha (t-\mathrm{epoch}),\end{array}$$(2)

where γ₁ and γ₂ are the velocity offsets of both instruments and α is the slope of the linear drift. The function δ₁(t) (respectively, δ₂(t)) is equal to one for the points taken by instrument 1 (respectively, instrument 2) and zero otherwise. The matrix φ_ℋ would thus be a n × 3 matrix whose columns would be the three explanatory time series φ_ℋ = (δ₁(t),δ₂(t),(t − epoch)) and the vector of parameters would be θ_ℋ = (γ₁, γ₂, α).

The enlarged model 𝒦(ν) is written as

$$\begin{array}{c}\hfill \mathcal{K}(\nu ):m_{\mathcal{K}}(\nu ,{\theta }_{\mathcal{K}})={\phi }_{\mathcal{K}}(\nu ){\theta }_{\mathcal{K}},\end{array}$$(3)

where θ_𝒦 = (θ_ℋ, θ) is the vector of size p + d of the parameters and φ_𝒦(ν) = (φ_H, φ(ν)) is a n × (p + d) matrix whose p first columns are those of φ_ℋ, and whose d last columns are functions of the frequency ν. Typically, d = 2, and the two additional columns are cos(νt) and sin(νt), but the theory developed by Davies (1977, 1987, 2002) and Baluev (2008) is more general.

We denote by ${\chi }_{\mathcal{H}}^2$ and ${\chi }_{\mathcal{K}}^2(\nu )$ the χ² of the residuals after a linear least squares fit with a covariance matrix C of the models ℋ and 𝒦(ν), respectively. Baluev (2008) assumed the noise to be independent (diagonal covariance matrix) and Gaussian and that the uncertainties of the measurements are known precisely (at least within a common factor). In this generalization, we do not assume the noise to be independent anymore, but we still assume the noise to be Gaussian with a known covariance matrix C (at least within a common factor). The covariance matrix C accounts for all sources of correlated and uncorrelated noise, such as intrinsic noise from the source or subsequent contamination by the Earth’s atmosphere or by the instrument.

In the general case, the periodogram is a function z(ν) = f(${\chi }_{\mathcal{H}}^2$, ${\chi }_{\mathcal{K}}^2(\nu )$). A general linear periodogram z is thus defined by the models ℋ and 𝒦(ν) and the function f. In the following, we consider the four definitions of the periodogram proposed by Baluev (2008):

$$\begin{array}{cc}\hfill z_0(\nu )=\frac{1}{2}({\chi }_{\mathcal{H}}^2-{\chi }_{\mathcal{K}}^2(\nu )),& z_1(\nu )=\frac{n_{\mathcal{H}}}{2}\frac{{\chi }_{\mathcal{H}}^2-{\chi }_{\mathcal{K}}^2(\nu )}{{\chi }_{\mathcal{H}}^2},\hfill \\ \hfill z_2(\nu )=\frac{n_{\mathcal{K}}}{2}\frac{{\chi }_{\mathcal{H}}^2-{\chi }_{\mathcal{K}}^2(\nu )}{{\chi }_{\mathcal{K}}^2(\nu )},& z_3(\nu )=\frac{n_{\mathcal{K}}}{2}ln\frac{{\chi }_{\mathcal{H}}^2}{{\chi }_{\mathcal{K}}^2(\nu )},\hfill \end{array}$$(4)

where n_ℋ = n − p and n_𝒦 = n − (p + d).

The widespread GLS periodogram (see Ferraz-Mello 1981; Zechmeister & Kürster 2009) is very close to the definition z₁ of the periodogram. Indeed, we have

$$\begin{array}{c}\hfill z_{\mathrm{GLS}}=\frac{{\chi }_{\mathcal{H}}^2-{\chi }_{\mathcal{K}}^2(\nu )}{{\chi }_{\mathcal{H}}^2}=\frac{2}{n_{\mathcal{H}}}{z}_1(\nu ),\end{array}$$(5)

and all the results obtained for z₁ are also valid for the GLS.

Once the periodogram is computed, it is useful to compute the p-value of the highest peak, or FAP, defined as Pr{max_νz(ν) ≥ Z | ℋ}, where Z is the value of the maximum peak of the periodogram computed on the data.

2.2. False alarm probability for periodograms with correlated noise

In this section, we provide analytical approximations of the FAP for the definitions of the periodogram of Eq. (4). Their precise derivation is provided in Appendix A.

The model ℋ is defined as in Eq. (1), where the n × p matrix φ_ℋ is user defined; for instance, it might include offsets and drifts. The model 𝒦 (Eq. (3)) is the horizontal concatenation of φ_ℋ and the two column vectors cos νt and sinνt (φ_𝒦(ν) = (φ_ℋ, cos νt, sinνt)).

The periodogram is computed in the range of frequencies ]0, ν_max]. The FAP is approximated by (see Baluev 2008)

$$\begin{array}{c}\hfill {\mathrm{FAP}}_{max}(Z,{\nu }_{\max })\approx 1-(1-{\mathrm{FAP}}_{\text{single}}(Z)){\mathit{e}}^{-\tau (Z,{\nu }_{\max })},\end{array}$$(6)

where analytical expressions for FAP_single and τ(Z, ν_max) are given in Table 1. These expressions depend on the rescaled frequency bandwidth W defined as

$$\begin{array}{c}\hfill W=\frac{{\nu }_{\max }}{2\pi }{T}_{\mathrm{eff}},\end{array}$$(7)

where T_eff is the effective time series length, which we approximate by (see Appendix A.2)

$$\begin{array}{c}\hfill T_{\mathrm{eff}}\approx \sqrt{4\pi }\sqrt{\frac{⟨\mathrm{\Pi }\ast \text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}{⟨\text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}-{\left(\frac{⟨\mathrm{\Sigma }\ast \text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}{2⟨\text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}\right)}^2}·\end{array}$$(8)

The n × n matrices Σ, Δ, and Π are defined as

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{i,j}& =t_i+t_j,\hfill \\ \hfill {\mathrm{\Delta }}_{i,j}& =t_i-t_j,\hfill \\ \hfill {\mathrm{\Pi }}_{i,j}& =t_i{t}_j,\hfill \end{array}$$(9)

and for two n × n matrices X and Y, X*Y is the Hadamard (or element-wise) product

$$\begin{array}{c}\hfill {(X\ast Y)}_{i,j}=X_{i,j}{Y}_{i,j},\end{array}$$(10)

and ⟨X⟩ is defined as

$$\begin{array}{c}\hfill ⟨X⟩={\displaystyle \sum _{i,j}}{C}_{i,j}^{-1}{X}_{i,j}.\end{array}$$(11)

The expression of the effective time series length found by Baluev (2008) in the white noise case can be derived from Eq. (8). Indeed in this case (diagonal covariance matrix C), Eq. (8) is simplified as

$$\begin{array}{c}\hfill T_{\mathrm{eff}}\approx \sqrt{4\pi (\overline{t^2}-{\overline{t}}^2)},\end{array}$$(12)

where $\overline{t}$, and $\overline{t^2}$ are weighted means with weights $C_{i,i}^{-1}/{\sum }_j{C}_{j,j}^{-1}$.

In Appendix A.2, we additionally provide approximations of T_eff in the low and high frequency limit, that is, ν_maxΔ ≪ 1 (Eq. (A.24)) and ν_maxΔ ≫ 1 (Eq. (A.25)).

2.3. Comparison of analytical FAP with Monte Carlo simulations

We validated our analytical estimate of the FAP (Eqs. (6)–(8), Table 1) by comparing it with Monte Carlo simulations. The Monte Carlo simulations are performed by generating a large set of random time series following the same distribution (same covariance matrix). We used the times of observation and error bars of the HARPS radial velocities of HD 136352 (Udry et al. 2019) to obtain a realistic temporal sampling and realistic covariance matrices. The HARPS radial velocities of HD 136352 consist of 648 points taken over almost 11 years (2004–2015) and spread over 238 distinct nights (about 2.7 points per night).

Our method is general and does not require a particular shape for the covariance matrix. However, for illustration purposes, we assume the covariance matrix to follow

$$\begin{array}{c}\hfill C_{i,j}={\delta }_{i,j}({\sigma }_i^2+{\sigma }_{\mathrm{jit}.}^2)+{\sigma }_{\exp .}^2{\mathit{e}}^{-|t_i-t_j|/{\tau }_{\exp .}}\end{array}$$(13)

and vary the values of the parameters (σ_jit., σ_exp., τ_exp.) to define four different noise models:

1. obs. (white noise): a diagonal matrix with observational error bars;

2. jit. (white noise): a diagonal matrix with observational error bars plus a jitter of 1 m s⁻¹;

3. daily exp. (correlated noise): observational error bars on the diagonal, plus an exponential decay of 1 m s⁻¹ with a timescale of 1 d;

4. monthly exp. (correlated noise): the same as daily exp. but with a timescale of 30 d.

The values of the parameters (σ_jit., σ_exp., τ_exp.) used for each noise model are summarized in Table 2. In the context of radial velocity time series, jitter terms might model both intrinsic noise from the star and instrumental noise, while exponential decay terms are often used to account for stellar noise (e.g., granulation and oscillation).

For a given covariance matrix C, we generated a synthetic radial velocity time series by randomly sampling from a normal distribution with covariance matrix C. We generate 10⁶ such random time series and compute a periodogram (with the correct covariance matrix) for each time series. The periodograms are computed in the range $]0,\frac{2\pi }{P_{\min }}]$ where P_min = 0.9 d, and with an instrumental offset γ adjusted for each frequency (p = 1, n_ℋ = n − 1, n_𝒦 = n − 3). Then, the distribution of the maximum of these periodograms allows us to estimate numerically the FAP.

The comparison between the analytical and numerical FAP is shown in Fig. 1. As explained by Baluev (2008), the analytical formula of the FAP is an upper bound that asymptotically (for low FAP levels) converges to the exact FAP. This is indeed what we observe in Fig. 1. For all the covariance matrices, the analytical and numerical estimates agree very well for FAP ≲ 0.1, and the analytical formula overestimates the FAP for FAP ≳ 0.1

Fig. 1. Comparison between analytical and numerical estimations of the FAP for four types of covariance matrix (see Sect. 2.3) and using the HARPS time series of HD 136352. The periodogram power is computed following the definition z1 of Eq. (4). The expectation of z1 is one (see Eq. (B.10)).

Since we used 10⁶ samples for the numerical estimation of the FAP, we could not reliably explore FAP levels below 10⁻⁴ owing to small numbers statistics. However, the analytical and numerical estimates agree very well down to 10⁻⁴, and the analytical approximation is expected to be even more accurate for lower FAP levels.

To sum up, the analytical estimate provides a very good approximation of the FAP in the range of most interest (FAP ≲ 0.1), and is conservative for higher FAP levels. Moreover, this analytical estimate is much faster to compute than Monte Carlo simulations (or other numerical methods), especially for low FAP levels. These properties make it very convenient to use in practical applications.

3. Sensitivity of periodogram to noise model

The FAP formula obtained in Sect. 2.2 provides an efficient and robust way to assess the significance of periodogram peaks when the covariance matrix is known. However, this is not the case in general, and we can often only make educated guesses about the shape of the covariance matrix. It is therefore necessary to explore the sensitivity of the periodogram to the noise model. To do so, we computed the expectation of the periodogram computed with an incorrect noise model. The analytical derivation of the periodogram expectation is described in Appendix B. In this section, we illustrate its use by exploring the effect of an incorrect noise model on the periodogram and its associated FAP.

We assume that the actual covariance matrix of the noise is C, while the periodogram is computed with an incorrect covariance matrix V. We used the same dataset (HARPS radial velocities of HD 136352) and the same noise models as in Sect. 2.3. We first chose a noise model among the four models of Sect. 2.3 (obs., jit., daily exp., and monthly exp.), which we considered the correct noise model (covariance matrix C). We then chose the incorrect noise model (covariance matrix V) among the three other models. The study is done with the definition z₁ of the periodogram (see Eq. (4)). We first computed the periodogram expectation following the analytical expression of Eq. (B.9) and we then estimated the FAP corresponding to the highest peak of this mean periodogram using the analytical formula of Eq. (6) with the incorrect covariance matrix (V).

The results are shown in Fig. 2 for all pairs of noise models. We observe in Fig. 2 that adding a jitter term (jit.) does not affect the results much compared to the model with the observational error bars alone (obs.), and vice versa. This is not surprising since all the data points were taken with the same instruments (HARPS), and thus have very similar error bars. Therefore, adding a common jitter term to all error bars does not much affect the relative weight of each measurement. Moreover, the definition z₁ of the periodogram (Eq. (4)) is not sensitive to the multiplication of all error bars by a common factor.

Fig. 2. Periodogram expectation for HD 136352 in the case of a wrong noise model (see Eq. (B.9)). Rows correspond to the true noise model, while columns correspond to the assumed wrong model. The definition z1 of the periodogram power is used (see Eq. (4)). If the noise model was correct, the expectation of the periodogram power would be uniformly 1 (see Eq. (B.10)). The red vertical lines highlight 0.5 d, 1 d, and 1 yr. For each periodogram, we provide the analytical estimate of the FAP (using the wrong noise model) corresponding to the highest peak of the periodogram expectation.

Using a correlated noise model (daily or monthly exp.) while the true noise is uncorrelated (obs. or jit.) remains conservative on the whole range of frequencies (𝔼(z₁(ν)) ≲ 1 for all ν). However, the periodogram level is very low at long periods, which means that the detection capability at long periods is strongly reduced. On the other hand, using an uncorrelated noise model (obs. or jit.) while the true noise model is correlated (daily or monthly exp.) can lead to spurious detections with very low FAP levels (down to 3.22 × 10⁻⁶, see Fig. 2). Underestimating the correlation timescale (using daily instead of monthly exp.) has a similar (but weaker) effect as using an uncorrelated model instead of a correlated model. Finally overestimating the correlation timescale (using monthly instead of daily exp.) reduces the capability to detect long periods, and might lead to spurious detections of short periods. In the two latter cases (underestimation and overestimation of the correlation timescale), the FAP remains very high (close to 1, i.e., non-significant detection).

Overall these results are not surprising but illustrate the possibility to investigate the sensitivity of the periodogram and its FAP with respect to the noise model using the formula for the periodogram expectation (Eq. (B.9)).

4. Conclusions

We present a generalization of the analytical estimate of Baluev (2008) (which was restricted to the white noise case) to the correlated noise case (see Sect. 2.2). We show that the “suggestive generalization” of Baluev (2013c) is valid in the low frequency limit (see Eq. (A.24)), but we find a more general expression (Eq. (8)) that is valid for all frequencies. We validate our analytical estimate against Monte Carlo simulations (see Sect. 2.3) and show that this analytical criterion is very efficient and accurate, provided that the covariance matrix of the noise is known (at least within a common factor).

In most cases, however, astronomical time series are contaminated by sources of correlated noise that are difficult to characterize fully, which results in an approximate modeling of the covariance matrix. We illustrate the sensitivity of the periodogram to the noise model, by deriving the expectation of a periodogram computed with an incorrect covariance matrix (see Sect. 3). This method allows us to visualize which parts of the periodogram are the most affected by a change in the noise model. For instance, we observe, as expected, that overestimating the correlation timescale of the noise tends to reduce the detection capability at long periods strongly, while underestimating this timescale can lead to spurious detections. Another way to visualize the sensitivity of the periodogram to the noise model would be to compute several periodograms on the same data with various noise models. Both approaches are complementary to better understand the features observed in the periodograms.

Several methods can be used to obtain a more realistic covariance matrix. First, a likelihood maximization can be performed to adjust some noise model parameters. This maximization can be performed once, before computing the least squares periodogram (with a fixed noise model). It can also be performed for each frequency by computing a likelihood periodogram instead of a least squares periodogram. In the white noise case, Baluev (2009) proposed a FAP estimate for a likelihood periodogram with a free error term (jitter) added in quadrature to the nominal error bars and adjusted for at each frequency. A generalization to the correlated noise case could probably also be achieved for this likelihood periodogram, but is beyond the scope of this article. Finally, a Bayesian approach could be used to compare different models (signal + noise), but with a much higher computational cost.

Acknowledgments

We thank the anonymous referee for his/her useful comments. We acknowledge financial support from the Swiss National Science Foundation (SNSF). This work has, in part, been carried out within the framework of the National Centre for Competence in Research PlanetS supported by SNSF.

References

Appendix A: Computation of the FAP in the correlated noise case

In this appendix, we extend the method of Baluev (2008) to obtain analytical FAP estimates in the correlated noise case. The main idea allowing the analytical approximation of the FAP with correlated noise is to perform a change of coordinates that yields independent Gaussian random variables. Then, the method described by Baluev (2008) can be applied on these new variables. The change of variables is described in Appendix A.1 and the derivation of the FAP estimate in Appendix A.2.

A.1. Change of random variables

Let us assume that the covariance matrix C of the noise is known (at least within a common factor). Then, under the null hypothesis (i.e., assuming the base model is correct), the time series is written as

$$\begin{array}{c}\hfill y={\phi }_{\mathcal{H}}{\theta }_{\mathcal{H},0}+ϵ,\end{array}$$(A.1)

where θ_ℋ, 0 is the true value of the parameters and the noise ϵ is Gaussian with zero mean and covariance C.

For a linear model φ_m and parameters θ_m (m = ℋ or 𝒦), the χ² is defined as

$$\begin{array}{cc}\hfill {\chi }^2({\theta }_m)& ={(y-{\phi }_m{\theta }_m)}^{\mathrm{T}}{C}^{-1}(y-{\phi }_m{\theta }_m)\hfill \\ \hfill & ={({\phi }_{\mathcal{H}}{\theta }_{\mathcal{H},0}-{\phi }_m{\theta }_m+ϵ)}^{\mathrm{T}}{C}^{-1}({\phi }_{\mathcal{H}}{\theta }_{\mathcal{H},0}-{\phi }_m{\theta }_m+ϵ).\hfill \end{array}$$(A.2)

The least squares estimate of the parameters is written as

$$\begin{array}{cc}\hfill {\widehat{\theta }}_m& ={\left({\phi }_m^{\mathrm{T}}{C}^{-1}{\phi }_m\right)}^{-1}{\phi }_m^{\mathrm{T}}{C}^{-1}y\hfill \\ \hfill & ={\theta }_{\mathcal{H},0}+{\left({\phi }_m^{\mathrm{T}}{C}^{-1}{\phi }_m\right)}^{-1}{\phi }_m^{\mathrm{T}}{C}^{-1}ϵ,\hfill \end{array}$$(A.3)

and the minimum χ² is thus

$$\begin{array}{cc}\hfill {\chi }_m^2& =\underset{{\theta }_m}{min}{\chi }^2({\theta }_m)={\chi }^2({\widehat{\theta }}_m)={(y-{\phi }_m\widehat{\theta })}^{\mathrm{T}}{C}^{-1}(y-{\phi }_m\widehat{\theta })\hfill \\ \hfill & =y^{\mathrm{T}}(C^{-1}-C^{-1}{\phi }_m{\left({\phi }_m^{\mathrm{T}}{C}^{-1}{\phi }_m\right)}^{-1}{\phi }_m^{\mathrm{T}}{C}^{-1})y\hfill \\ \hfill & ={ϵ}^{\mathrm{T}}(C^{-1}-C^{-1}{\phi }_m{\left({\phi }_m^{\mathrm{T}}{C}^{-1}{\phi }_m\right)}^{-1}{\phi }_m^{\mathrm{T}}{C}^{-1})ϵ.\hfill \end{array}$$(A.4)

Let us now perform the following change of coordinates:

$$\begin{array}{cc}& z=L^{-1}y,\hfill \\ \hfill & \eta =L^{-1}ϵ,\hfill \\ \hfill & {\psi }_m=L^{-1}{\phi }_m,\hfill \end{array}$$(A.5)

where C = LL^T is the Cholesky decomposition of the covariance matrix. Since we assumed ϵ to follow a Gaussian law with zero mean and covariance C, η follows a Gaussian law with zero mean and covariance 𝟙. The random variables η are thus independent Gaussian variables. In these new variables, the χ² is simply rewritten as

$$\begin{array}{c}\hfill {\chi }^2(\theta )={(z-{\psi }_m{\theta }_m)}^{\mathrm{T}}(z-{\psi }_m{\theta }_m),\end{array}$$(A.6)

the least squares estimate of the parameters is rewritten as

$$\begin{array}{c}\hfill {\widehat{\theta }}_m={\theta }_{\mathcal{H},0}+{\left({\psi }_m^{\mathrm{T}}{\psi }_m\right)}^{-1}{\psi }_m^{\mathrm{T}}\eta ,\end{array}$$(A.7)

and the minimum χ² as

$$\begin{array}{c}\hfill {\chi }_m^2={\eta }^{\mathrm{T}}(\mathbb{1}-{\psi }_m{\left({\psi }_m^{\mathrm{T}}{\psi }_m\right)}^{-1}{\psi }_m^{\mathrm{T}})\eta ,\end{array}$$(A.8)

which follows a χ² law with n_m degrees of freedom (n_ℋ = n − p, n_𝒦 = n − (p + d)). Therefore, the initial problem of analyzing a time series y with covariance matrix C, base model φ_ℋ, and enlarged models φ_𝒦 = (φ_H, φ(ν)) is equivalent to analyzing the time series z, with covariance matrix 𝟙, base model ψ_ℋ, and enlarged model ψ_𝒦 = (ψ_H, ψ(ν)). However, if φ(ν) was the sine and cosine at frequency ν, this is no longer the case for ψ(ν). Nevertheless, the theory of Baluev (2008) is very general, and does not assume a particular shape for this matrix, except for the final application to the least squares periodogram. We thus follow the method proposed by Baluev (2008), and only change the hypothesis on the shape of the enlarged model matrix.

A.2. Analytical FAP estimate

The FAP can be bounded by (see Baluev 2008, Eq. (5))

$$\begin{array}{c}\hfill {\text{FAP}}_{max}(Z,{\nu }_{\max })\le {\text{FAP}}_{\text{single}}(Z)+\tau (Z,{\nu }_{\max }),\end{array}$$(A.9)

and approximated by (see Baluev 2008, Eq. (6))

$$\begin{array}{c}\hfill {\text{FAP}}_{max}(Z,{\nu }_{\max })\approx 1-(1-{\text{FAP}}_{\text{single}}(Z)){\mathit{e}}^{-\tau (Z,{\nu }_{\max })},\end{array}$$(A.10)

where Z is the maximum periodogram power, FAP_single(Z) is the FAP in the case in which the frequency ν of the putative additional signal is fixed, τ(Z, ν_max) is the expectation of the number of up-crossings of the level Z by the periodogram (see Baluev 2008).

Computing FAP_single(Z) and τ(Z, ν_max) requires us to specify the definition of the periodogram z(ν). Baluev (2008) proposed several definitions and derived the corresponding formulas for FAP_single(Z) and τ(Z, ν_max). These results are summarized in Table 1 for d = 2. The general case (d not necessarily equal to 2 and other definitions of z(ν)) is provided in Baluev (2008), Appendix B.

For the definitions of the periodogram of Eq. (4) and assuming d = 2, the only quantity left to compute is the factor W, which is the rescaled frequency bandwidth, defined as (see Baluev 2008)

$$\begin{array}{c}\hfill W=\frac{A({\nu }_{\max })}{2{\pi }^{3/2}},\end{array}$$(A.11)

where

$$\begin{array}{cc}\hfill A({\nu }_{\max })=& {\displaystyle {\int }_0^{{\nu }_{\max }}}\mathrm{d}\nu {\displaystyle {\int }_{x^2<1}}\frac{\sqrt{x^{\mathrm{T}}M(\nu )x}}{x^{\mathrm{T}}x}\mathrm{d}x\hfill \\ \hfill & \le 2\pi {\displaystyle {\int }_0^{{\nu }_{\max }}}\sqrt{\frac{\text{tr}(M(\nu ))}{2}}\mathrm{d}\nu .\hfill \end{array}$$(A.12)

The 2 × 2 matrix M(ν) is defined as follows (with x′ = ∂x/∂ν):

$$\begin{array}{cc}& Q={\psi }^{\mathrm{T}}\psi ={\phi }^{\mathrm{T}}{C}^{-1}\phi ,S={\psi }^{\mathrm{T}}{\psi }^{\prime }={\phi }^{\mathrm{T}}{C}^{-1}{\phi }^{\prime },\hfill \\ \hfill & R={\psi }^{\prime \mathrm{T}}{\psi }^{\prime }={\phi }^{\prime \mathrm{T}}{C}^{-1}{\phi }^{\prime },\hfill \\ \hfill & Q_{\mathcal{H}}={\psi }_{\mathcal{H}}^{\mathrm{T}}\psi ={\phi }_{\mathcal{H}}^{\mathrm{T}}{C}^{-1}\phi ,S_{\mathcal{H}}={\psi }_{\mathcal{H}}^{\mathrm{T}}{\psi }^{\prime }={\phi }_{\mathcal{H}}^{\mathrm{T}}{C}^{-1}{\phi }^{\prime },\hfill \\ \hfill & Q_{\mathcal{H},\mathcal{H}}={\psi }_{\mathcal{H}}^{\mathrm{T}}{\psi }_{\mathcal{H}}={\phi }_{\mathcal{H}}^{\mathrm{T}}{C}^{-1}{\phi }_{\mathcal{H}},\hfill \\ \hfill & \stackrel{\sim }{Q}=Q-Q_{\mathcal{H}}^{\mathrm{T}}{Q}_{\mathcal{H},\mathcal{H}}^{-1}{Q}_{\mathcal{H}},\stackrel{\sim }{S}=S-Q_{\mathcal{H}}^{\mathrm{T}}{Q}_{\mathcal{H},\mathcal{H}}^{-1}{S}_{\mathcal{H}},\hfill \\ \hfill & \stackrel{\sim }{R}=R-S_{\mathcal{H}}^{\mathrm{T}}{Q}_{\mathcal{H},\mathcal{H}}^{-1}{S}_{\mathcal{H}},\hfill \\ \hfill & M={\stackrel{\sim }{Q}}^{-1}(\stackrel{\sim }{R}-{\stackrel{\sim }{S}}^{\mathrm{T}}{\stackrel{\sim }{Q}}^{-1}\stackrel{\sim }{S}).\hfill \end{array}$$(A.13)

Baluev (2008) also defined the effective time series length as

$$\begin{array}{c}\hfill T_{\mathrm{eff}}=\frac{A({\nu }_{\max })}{\sqrt{\pi }{\nu }_{\max }},\end{array}$$(A.14)

such that

$$\begin{array}{c}\hfill W=\frac{{\nu }_{\max }}{2\pi }{T}_{\mathrm{eff}}.\end{array}$$(A.15)

From Eqs. (A.12) and (A.14), we obtain

$$\begin{array}{c}\hfill T_{\mathrm{eff}}=\frac{1}{\sqrt{\pi }}\overline{{\displaystyle {\int }_{x^2<1}}\frac{\sqrt{x^{\mathrm{T}}M(\nu )x}}{x^{\mathrm{T}}x}\mathrm{d}x}\le \overline{\sqrt{2\pi \text{tr}(M(\nu ))}},\end{array}$$(A.16)

where $\overline{x}$ is the mean of x(ν) over the frequency range ]0, ν_max]. As noted by Baluev (2008), the inequality in Eqs. (A.12) and (A.16) is very sharp in practical applications. In particular, it saturates (i.e., becomes an equality) when the eigenvalues of M(ν) are equal. This expression can be evaluated numerically by sampling the frequency over the interval ]0, ν_max], and computing M(ν) according to Eq. (A.13) for each frequency ν. The cost of evaluating T_eff is of the same order of magnitude as computing the periodogram itself. It is therefore much more efficient than performing Monte Carlo simulations. However, this cost is not negligible compared to the periodogram cost and analytical approximations might be useful.

We now specify the expression of T_eff for φ = (cos(νt), sin(νt)). Replacing φ in the definitions of Q, S, and R, we find

$$\begin{array}{cc}\hfill Q& =\frac{1}{2}\left(\begin{array}{cc}⟨cos\nu \mathrm{\Sigma }+cos\nu \mathrm{\Delta }⟩& ⟨sin\nu \mathrm{\Sigma }⟩\\ ⟨sin\nu \mathrm{\Sigma }⟩& ⟨cos\nu \mathrm{\Delta }-cos\nu \mathrm{\Sigma }⟩\end{array}\right),\hfill \\ \hfill S& =\frac{1}{4}\left(\begin{array}{cc}-⟨\mathrm{\Sigma }\ast sin\nu \mathrm{\Sigma }+\mathrm{\Delta }\ast sin\nu \mathrm{\Delta }⟩& ⟨\mathrm{\Sigma }\ast (cos\nu \mathrm{\Sigma }+cos\nu \mathrm{\Delta })⟩\\ ⟨\mathrm{\Sigma }\ast (cos\nu \mathrm{\Sigma }-cos\nu \mathrm{\Delta })⟩& ⟨\mathrm{\Sigma }\ast sin\nu \mathrm{\Sigma }-\mathrm{\Delta }\ast sin\nu \mathrm{\Delta }⟩\end{array}\right),\hfill \\ \hfill R& =\frac{1}{2}\left(\begin{array}{cc}⟨\mathrm{\Pi }\ast (cos\nu \mathrm{\Delta }-cos\nu \mathrm{\Sigma })⟩& -⟨\mathrm{\Pi }\ast sin\nu \mathrm{\Sigma }⟩\\ -⟨\mathrm{\Pi }\ast sin\nu \mathrm{\Sigma }⟩& ⟨\mathrm{\Pi }\ast (cos\nu \mathrm{\Sigma }+cos\nu \mathrm{\Delta })⟩\end{array}\right),\hfill \end{array}$$(A.17)

where  *  denotes the Hadamard (or elementwise) product,

$$\begin{array}{cc}\hfill ⟨X⟩& ={\displaystyle \sum _{i,j}}{C}_{i,j}^{-1}{X}_{i,j},\hfill \\ \hfill {\mathrm{\Sigma }}_{i,j}& =t_i+t_j,\hfill \\ \hfill {\mathrm{\Delta }}_{i,j}& =t_i-t_j,\hfill \\ \hfill {\mathrm{\Pi }}_{i,j}& =t_i{t}_j.\hfill \end{array}$$(A.18)

We follow Baluev (2008) and neglect aliasing effects. In this approximation all the terms containing sine or cosine of νΣ average out. The terms containing sinνΔ can also be neglected. Indeed, in the low frequency limit (νΔ ≪ 1), the terms in sinνΔ vanish, while in the high frequency limit (νΔ ≫ 1), the terms in sinνΔ average out. We thus obtain

$$\begin{array}{cc}\hfill Q& \approx \frac{1}{2}⟨cos\nu \mathrm{\Delta }⟩\mathbb{1},\hfill \\ \hfill S& \approx \frac{1}{4}⟨\mathrm{\Sigma }\ast cos\nu \mathrm{\Delta }⟩J,\hfill \\ \hfill R& \approx \frac{1}{2}⟨\mathrm{\Pi }\ast cos\nu \mathrm{\Delta }⟩\mathbb{1},\hfill \end{array}$$(A.19)

where J is the antisymmetric matrix

$$\begin{array}{c}\hfill J=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right).\end{array}$$(A.20)

As in Baluev (2008), we additionally assume that ψ(ν) is orthogonal to ψ_ℋ for all ν. As a consequence, Q̃ = Q, S̃ = S, R̃ = R, and

$$\begin{array}{c}\hfill M(\nu )=Q^{-1}(R-S^{\mathrm{T}}{Q}^{-1}S).\end{array}$$(A.21)

Replacing Eq. (A.19) in Eq. (A.21), we find

$$\begin{array}{c}\hfill M(\nu )\approx (\frac{⟨\mathrm{\Pi }\ast cos\nu \mathrm{\Delta }⟩}{⟨cos\nu \mathrm{\Delta }⟩}-{\left(\frac{⟨\mathrm{\Sigma }\ast cos\nu \mathrm{\Delta }⟩}{2⟨cos\nu \mathrm{\Delta }⟩}\right)}^2)\mathbb{1}.\end{array}$$(A.22)

The two eigenvalues of M are thus equal in this approximation, and we can approximate T_eff with (see Eq. (A.16))

$$\begin{array}{c}\hfill T_{\mathrm{eff}}\approx \sqrt{4\pi }\overline{\sqrt{\frac{⟨\mathrm{\Pi }\ast cos\nu \mathrm{\Delta }⟩}{⟨cos\nu \mathrm{\Delta }⟩}-{\left(\frac{⟨\mathrm{\Sigma }\ast cos\nu \mathrm{\Delta }⟩}{2⟨cos\nu \mathrm{\Delta }⟩}\right)}^2}}·\end{array}$$(A.23)

In the low frequency limit ν_maxΔ ≪ 1, the cosines can all be replaced by 1, and we obtain

$$\begin{array}{cc}\hfill T_{}\mathrm{eff},low& \approx \sqrt{4\pi }\sqrt{\frac{⟨\mathrm{\Pi }⟩}{⟨1⟩}-{\left(\frac{⟨\mathrm{\Sigma }⟩}{2⟨1⟩}\right)}^2}\hfill \\ \hfill & \approx \sqrt{4\pi }\sqrt{\frac{t^{\mathrm{T}}{C}^{-1}t}{u^{\mathrm{T}}{C}^{-1}u}-{\left(\frac{u^{\mathrm{T}}{C}^{-1}t}{u^{\mathrm{T}}{C}^{-1}u}\right)}^2},\hfill \end{array}$$(A.24)

where u is the vector of size n filled with ones. This expression was proposed by Baluev (2013c) as a “suggestive generalization” of the results of Baluev (2008) to correlated noise. However, Baluev (2013c) advocates against its use since it was not proved rigorously. Moreover, this expression is not valid in the general case but only in the low frequency limit.

In the high frequency limit ν_maxΔ ≫ 1, the cosines average out in Eq. (A.23), except on the diagonal (where they are equal to 1), and we find

$$\begin{array}{c}\hfill T_{}\mathrm{eff},high\approx \sqrt{4\pi }\sqrt{\frac{w^{\mathrm{T}}(t\ast t)}{w^{\mathrm{T}}u}-{\left(\frac{w^{\mathrm{T}}t}{w^{\mathrm{T}}u}\right)}^2},\end{array}$$(A.25)

where w = diag(C⁻¹).

Finally, for any frequency ν_max, we can approximate (at first order) the integral of Eq. (A.23) by replacing each sum ⟨X(ν)⟩ by its average $\overline{⟨X⟩}$ over the range ]0, ν_max]. We have

$$\begin{array}{cc}& \overline{⟨cos\nu \mathrm{\Delta }⟩}=⟨\text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩,\hfill \\ \hfill & \overline{⟨\mathrm{\Sigma }\ast cos\nu \mathrm{\Delta }⟩}=⟨\mathrm{\Sigma }\ast \text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩,\hfill \\ \hfill & \overline{⟨\mathrm{\Pi }\ast cos\nu \mathrm{\Delta }⟩}=⟨\mathrm{\Pi }\ast \text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩,\hfill \end{array}$$(A.26)

and thus

$$\begin{array}{c}\hfill T_{\mathrm{eff}}\approx \sqrt{4\pi }\sqrt{\frac{⟨\mathrm{\Pi }\ast \text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}{⟨\text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}-{\left(\frac{⟨\mathrm{\Sigma }\ast \text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}{2⟨\text{sinc}{\nu }_{\max }\mathrm{\Delta }⟩}\right)}^2}·\end{array}$$(A.27)

Equations (A.24) and (A.25) can also be derived directly from Eq. (A.27) in the low frequency and high frequency approximations. Moreover, as explained in Sect. 2.2, the expression found by Baluev (2008) in the white noise case can also be derived from Eq. (A.27) by using the fact that the covariance matrix is diagonal.

In practical applications, all estimations of T_eff (numerical evaluation of Eqs. (A.16) or (A.23), or directly using Eqs. (A.24), (A.25), and (A.27)) yield similar results. Moreover, as noted by Baluev (2008) in the case of independent Gaussian noise, T_eff is often of the same order of magnitude as the total time span of the time series (T_span = max(t)−min(t)).

Appendix B: Periodogram expectation

In this appendix, we show how to obtain an analytical estimate of the expectation of a periodogram computed with an incorrect noise model. We assume that the actual covariance matrix C of the noise is not known and that the χ² and periodograms are computed using an incorrect covariance matrix V. Under the null hypothesis (model ℋ), the time series still follows Eq. (A.1) but the χ² becomes

$$\begin{array}{cc}\hfill {\chi }^2({\theta }_m)& ={(y-{\phi }_m{\theta }_m)}^{\mathrm{T}}{V}^{-1}(y-{\phi }_m{\theta }_m)\hfill \\ \hfill & ={({\phi }_{\mathcal{H}}{\theta }_{\mathcal{H},0}-{\phi }_m{\theta }_m+ϵ)}^{\mathrm{T}}{V}^{-1}({\phi }_{\mathcal{H}}{\theta }_{\mathcal{H},0}-{\phi }_m{\theta }_m+ϵ).\hfill \end{array}$$(B.1)

The least squares estimate of the parameters becomes

$$\begin{array}{c}\hfill {\widehat{\theta }}_m={\theta }_{\mathcal{H},0}+{\left({\phi }_m^{\mathrm{T}}{V}^{-1}{\phi }_m\right)}^{-1}{\phi }_m^{\mathrm{T}}{V}^{-1}ϵ,\end{array}$$(B.2)

and the minimum χ² is thus

$$\begin{array}{c}\hfill {\chi }_m^2={ϵ}^{\mathrm{T}}(V^{-1}-V^{-1}{\phi }_m{\left({\phi }_m^{\mathrm{T}}{V}^{-1}{\phi }_m\right)}^{-1}{\phi }_m^{\mathrm{T}}{V}^{-1})ϵ.\end{array}$$(B.3)

We introduce a change of coordinates that is slightly different from Appendix A.1 (Eq. (A.5)), i.e.,

$$\begin{array}{cc}& \eta =L^{-1}ϵ,\hfill \\ \hfill & {\zeta }_m=M^{-1}{\phi }_m,\hfill \\ \hfill & N=M^{-1}L,\hfill \end{array}$$(B.4)

where C = LL^T and V = MM^T are the Cholesky decompositions of the true and assumed covariance matrices, and η is a vector of independent, centered, and reduced Gaussian random variables. In these coordinates, the minimum χ² (Eq. (B.3)) is rewritten as

$$\begin{array}{c}\hfill {\chi }_m^2={\eta }^{\mathrm{T}}{N}^{\mathrm{T}}(\mathbb{1}-P_{V,m})N\eta ,\end{array}$$(B.5)

where

$$\begin{array}{c}\hfill P_{V,m}={\zeta }_m{\left({\zeta }_m^{\mathrm{T}}{\zeta }_m\right)}^{-1}{\zeta }_m^{\mathrm{T}}\end{array}$$(B.6)

is the projection matrix on the subspace of ℝⁿ defined by the vectors of ζ_m. The expectation of the minimum χ² with the wrong covariance matrix V is thus

$$\begin{array}{cc}\hfill {\mu }_m& =\mathbb{E}({\chi }_m^2)=\text{tr}\left(N^{\mathrm{T}}(\mathbb{1}-P_{V,m})N\right)\hfill \\ \hfill & =\text{tr}\left((\mathbb{1}-P_{V,m})C_V\right),\hfill \end{array}$$(B.7)

where C_V = NN^T = M⁻¹CM^-T. In the case V = C, we have C_V = 𝟙, and we deduce

$$\begin{array}{cc}\hfill {\mu }_{\mathcal{H}}& =n_{\mathcal{H}},\hfill \\ \hfill {\mu }_{\mathcal{K}}& =n_{\mathcal{K}}.\hfill \end{array}$$(B.8)

B.1. First order formula

At first order, the expectation of the periodogram can be obtained by replacing ${\chi }_{\mathcal{H}}^2$ (respectively, ${\chi }_{\mathcal{K}}^2$) by its expectation μ_ℋ (respectively, μ_𝒦) in the definition of the periodogram (Eq. (4)). We find

$$\begin{array}{cc}\hfill \mathbb{E}(z_0(\nu ))=\frac{1}{2}({\mu }_{\mathcal{H}}-{\mu }_{\mathcal{K}}(\nu )),& \mathbb{E}(z_1(\nu ))\approx \frac{n_{\mathcal{H}}}{2}\frac{{\mu }_{\mathcal{H}}-{\mu }_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{H}}},\hfill \\ \hfill \mathbb{E}(z_2(\nu ))\approx \frac{n_{\mathcal{K}}}{2}\frac{{\mu }_{\mathcal{H}}-{\mu }_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{K}}(\nu )},& \mathbb{E}(z_3(\nu ))\approx \frac{n_{\mathcal{K}}}{2}ln\frac{{\mu }_{\mathcal{H}}}{{\mu }_{\mathcal{K}}(\nu )}·\hfill \end{array}$$(B.9)

In the case V = C (the actual covariance matrix is known), we deduce

$$\begin{array}{c}\hfill \mathbb{E}(z_i(\nu ))\approx \frac{d}{2},\end{array}$$(B.10)

for i = 0, …, 3, and for all frequencies ν. However, in the case V ≠ C (wrong noise model), the periodogram expectation can significantly depart from d/2 and depends on the frequency ν.

B.2. Higher order formulas

Higher order estimates can also be obtained by developing Eq. (4) in power series of ${\chi }_m^2-{\mu }_m$ and computing higher order momenta of ${\chi }_{\mathcal{H}}^2$, ${\chi }_{\mathcal{K}}^2$. We provide more details in the following, however, the first order formula already yields very accurate results, and we thus adopt it in our study.

We introduce the random variables $X_m={\chi }_m^2-{\mu }_m$ (m = ℋ,  𝒦), which we assume to be small with respect to μ_m. We then develop the periodogram power of Eq. (4) in power series of X_m. For instance, at second order we obtain

$$\begin{array}{cc}\hfill z_0(\nu )=& \frac{1}{2}({\mu }_{\mathcal{H}}-{\mu }_{\mathcal{K}}(\nu )+X_{\mathcal{H}}-X_{\mathcal{K}}(\nu )),\hfill \\ \hfill z_1(\nu )=& \frac{n_{\mathcal{H}}}{2}(1-\frac{{\mu }_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{H}}}(1+\frac{X_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{K}}(\nu )}-\frac{X_{\mathcal{H}}}{{\mu }_{\mathcal{H}}}\hfill \\ \hfill & -\frac{X_{\mathcal{H}}{X}_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{H}}{\mu }_{\mathcal{K}}(\nu )}+{\left(\frac{X_{\mathcal{H}}}{{\mu }_{\mathcal{H}}}\right)}^2))+\mathcal{O}\left(X^3\right),\hfill \\ \hfill z_2(\nu )=& \frac{n_{\mathcal{K}}}{2}(-1+\frac{{\mu }_{\mathcal{H}}}{{\mu }_{\mathcal{K}}(\nu )}(1+\frac{X_{\mathcal{H}}}{{\mu }_{\mathcal{H}}}-\frac{X_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{K}}(\nu )}\hfill \\ \hfill & -\frac{X_{\mathcal{H}}{X}_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{H}}{\mu }_{\mathcal{K}}(\nu )}+{\left(\frac{X_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{K}}(\nu )}\right)}^2))+\mathcal{O}\left(X^3\right),\hfill \\ \hfill z_3(\nu )=& \frac{n_{\mathcal{K}}}{2}(ln\frac{{\mu }_{\mathcal{H}}}{{\mu }_{\mathcal{K}}(\nu )}+\frac{X_{\mathcal{H}}}{{\mu }_{\mathcal{H}}}-\frac{X_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{K}}(\nu )}\hfill \\ \hfill & -\frac{1}{2}{\left(\frac{X_{\mathcal{H}}}{{\mu }_{\mathcal{H}}}\right)}^2+\frac{1}{2}{\left(\frac{X_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{K}}(\nu )}\right)}^2)+\mathcal{O}\left(X^3\right).\hfill \end{array}$$(B.11)

The expectation of the periodogram is thus (at second order)

$$\begin{array}{cc}\hfill \mathbb{E}(z_0(\nu ))& =\frac{1}{2}({\mu }_{\mathcal{H}}-{\mu }_{\mathcal{K}}(\nu )),\hfill \\ \hfill \mathbb{E}(z_1(\nu ))& \approx \frac{n_{\mathcal{H}}}{2}(1-\frac{{\mu }_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{H}}}+\frac{\text{cov}({\chi }_{\mathcal{H}}^2,{\chi }_{\mathcal{K}}^2(\nu ))}{{\mu }_{\mathcal{H}}^2}-\frac{\text{var}({\chi }_{\mathcal{H}}^2){\mu }_{\mathcal{K}}(\nu )}{{\mu }_{\mathcal{H}}^3}),\hfill \\ \hfill \mathbb{E}(z_2(\nu ))& \approx \frac{n_{\mathcal{K}}}{2}(\frac{{\mu }_{\mathcal{H}}}{{\mu }_{\mathcal{K}}(\nu )}-1+\frac{\text{var}({\chi }_{\mathcal{K}}^2(\nu )){\mu }_{\mathcal{H}}}{{\mu }_{\mathcal{K}}{(\nu )}^3}-\frac{\text{cov}({\chi }_{\mathcal{H}}^2,{\chi }_{\mathcal{K}}^2(\nu ))}{{\mu }_{\mathcal{K}}{(\nu )}^2}),\hfill \\ \hfill \mathbb{E}(z_3(\nu ))& \approx \frac{n_{\mathcal{K}}}{2}(ln\frac{{\mu }_{\mathcal{H}}}{{\mu }_{\mathcal{K}}(\nu )}+\frac{\text{var}({\chi }_{\mathcal{K}}^2(\nu ))}{2{\mu }_{\mathcal{K}}{(\nu )}^2}-\frac{\text{var}({\chi }_{\mathcal{H}}^2)}{2{\mu }_{\mathcal{H}}^2}),\hfill \end{array}$$(B.12)

where μ_ℋ, μ_𝒦(ν) are computed according to Eq. (B.7), and

$$\begin{array}{cc}\hfill \text{cov}({\chi }_m^2,{\chi }_{m^{\prime }}^2)& =\text{tr}\left(N^{\mathrm{T}}(\mathbb{1}-P_{V,m})N{N}^{\mathrm{T}}(\mathbb{1}-P_{V,m^{\prime }})N\right)\hfill \\ \hfill & =\text{tr}\left((\mathbb{1}-P_{V,m})C_V(\mathbb{1}-P_{V,m^{\prime }})C_V\right).\hfill \end{array}$$(B.13)

All Tables

All Figures

	Fig. 1. Comparison between analytical and numerical estimations of the FAP for four types of covariance matrix (see Sect. 2.3) and using the HARPS time series of HD 136352. The periodogram power is computed following the definition z1 of Eq. (4). The expectation of z1 is one (see Eq. (B.10)).
In the text

Fig. 2. Periodogram expectation for HD 136352 in the case of a wrong noise model (see Eq. (B.9)). Rows correspond to the true noise model, while columns correspond to the assumed wrong model. The definition z1 of the periodogram power is used (see Eq. (4)). If the noise model was correct, the expectation of the periodogram power would be uniformly 1 (see Eq. (B.10)). The red vertical lines highlight 0.5 d, 1 d, and 1 yr. For each periodogram, we provide the analytical estimate of the FAP (using the wrong noise model) corresponding to the highest peak of the periodogram expectation.

In the text