Bias and robustness of eccentricity estimates from radial velocity data

TL;DR

This paper investigates the reliability of eccentricity estimates from radial velocity data, emphasizing the importance of proper noise modeling and convergence checks, and introduces methods to distinguish true planetary signals from artifacts.

Contribution

It provides a comprehensive analysis of eccentricity estimation robustness, proposing criteria for trustworthy results and methods to detect model inaccuracies in radial velocity data.

Findings

Reliable eccentricity estimates require proper noise modeling and convergence verification.

Strong peaks in the noise power spectrum can mimic eccentric signals.

Bayes factors help distinguish true planetary eccentricities from artifacts.

Abstract

Eccentricity is a parameter of particular interest as it is an informative indicator of the past of planetary systems. It is however not always clear whether the eccentricity fitted on radial velocity data is real or if it is an artefact of an inappropriate modelling. In this work, we address this question in two steps: we first assume that the model used for inference is correct and present interesting features of classical estimators. Secondly, we study whether the eccentricity estimates are to be trusted when the data contain incorrectly modelled signals, such as missed planetary companions, non Gaussian noises, correlated noises with unknown covariance, etc. Our main conclusion is that data analysis via posterior distributions, with a model including a free error term gives reliable results provided two conditions. First, convergence of the numerical methods needs to be ascertained. Secondly, the noise power spectrum should not have a particularly strong peak at the semi period of the planet of interest. As a consequence, it is difficult to determine if the signal of an apparently eccentric planet might be due to another inner companion in 2:1 mean motion resonance. We study the use of Bayes factors to disentangle these cases. Finally, we suggest methods to check if there are hints of an incorrect model in the residuals. We show on simulated data the performance of our methods and comment on the eccentricities of Proxima b and 55 Cnc f.

Tables

Table 1. Table 1: List of symbols.

$𝜽$	Vector of parameters
${𝜽}_t$	True value of the vector of parameters
$\widehat{𝜽}$	Estimator of $𝜽$
$\text{f}(\text{t},𝜽)$	Deterministic model sampled at times $\text{t}={\left(t_k\right)}_{k=1..N}$ and of parameters $𝜽$ (orbital parameters plus possibly offset, trend…)
$p\left(𝜽\right)$	Probability density of $𝜽$, $\mathrm{Pr}\left\{𝜽\in \Theta \right\}={\int }_{\Theta }p\left(𝜽\right)𝑑𝜽$ for some measurable set $\Theta$
$𝔼\left\{𝜽\right\}$	Mathematical expectancy of the random variable $𝜽$
$k,h$	$k=e\cos \omega$, $h=e\sin \omega$

Table 2. Table 2: Dimension of the models with fixed eccentricity as a function of the eccentricity, GJ876 measurement times

Eccentricity	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	0.999
Dimension of ${ℳ}_{e,P}$	3	4	4	6	8	10	14	24	46	91

Table 3. Table 3: Root mean squared error ( MSE MSE \sqrt{\mathrm{MSE}} ) and mean absolute error (MAE) averaged over all types of noises of several estimators: maximum likelihood ( e ML subscript 𝑒 ML e_{\mathrm{ML}} , max. a posteriori e MAP subscript 𝑒 MAP e_{\mathrm{MAP}} , mode, mean and median of the eccentricity posterior e max subscript 𝑒 max e_{\mathrm{max}} , e mean subscript 𝑒 mean e_{\mathrm{mean}} and e med subscript 𝑒 med e_{\mathrm{med}} for the K = 3.5 𝐾 3.5 K=3.5 m.s -1 simulation.

Estimator	$e_{\mathrm{ML}}$	$e_{\mathrm{MAP}}$	$e_{\max }$	$e_{\mathrm{mean}}$	$e_{\mathrm{med}}$
$\sqrt{\mathrm{MSE}}$	0.1384	0.1362	0.1613	0.1088	0.1117
MAE	0.1071	0.1062	0.1239	0.0842	0.0849

Table 4. Table 4: Mean absolute and mean square error (MAE and MSE) of the estimate when the eccentricities are generated with a Beta distributions, and the analysis is done with the same Beta distribution as prior (correct model) or done with a uniform prior (incorrect model)

Estimator	prior	K=2 m.s-1	K=3.5 m.s-1	K= 5 m.s-1
$\sqrt{\mathrm{MSE}}$	correct	0.1325	0.0994	0.0791
$\sqrt{\mathrm{MSE}}$	incorrect	0.1530	0.1118	0.0899
MAE	correct	0.1073	0.0803	0.0601
MAE	incorrect	0.1188	0.0851	0.0713

Table 5. Table 5: Priors used for the numerical experiments.

Parameter	prior
$K$	Uniform on $[0,10000]$
$P$	$1/P$ uniform on $[0,20]$
$e$	Beta(a=0.867,b=3.03) as in Kipping (2014)
$\omega$	uniform on $[0,2\pi ]$
$M_0$	uniform on $[0,2\pi ]$

Table 6. Table 6: Priors used for the Bayes factor of eccentric model vs 2:1 mean motion resonance. The symbol T obs subscript 𝑇 obs T_{\mathrm{obs}} denotes the observation time, P 0 = 60 subscript 𝑃 0 60 P_{0}=60 days and α = 0.1 𝛼 0.1 \alpha=0.1 .

Parameter	prior
${\sigma }_J^2$	Uniform on $[0,100]$ m/s
offset	Uniform on $[-100,100]$
$K,K_1$	Uniform in $\ln K$ on $[-1,9]$
$P,P_1$	$1/P$ uniform on $[1/P_0-1/T_{\mathrm{obs}},1/P_0+1/T_{\mathrm{obs}}]$
$e$	Uniform on $[0,1]$
$\omega$	uniform on $[0,2\pi ]$
$M_0$	uniform on $[0,2\pi ]$
$K_2$	Uniform in $\ln K$ on $[-1,9]$
$P_2$	$1/P$ uniform on $[\left(1-\alpha \right)\frac{2}{P_0},\left(1+\alpha \right)\frac{2}{P_0}]$
${\varphi }_1$	uniform on $[0,2\pi ]$
${\varphi }_2$	uniform on $[0,2\pi ]$

Table 7. Table 7: Error, mean squared error and noise level estimate of the eccentricity estimates for Gaussian and non Gaussian noises .

Noise type	Gauss- ian	Stu- dent $T$ 3	Stu- dent $T$ 4	Uni- form	Expo- nential	Poi- sson
Mean error	0.067	0.066	0.0672	0.067	0.066	0.067
$\sqrt{\mathrm{MSE}}$	0.084	0.094	0.090	0.085	0.083	0.084
Noise level	0.951	0.904	0.938	0.952	0.940	0.947

Equations

$$\text{y}(\text{t})=\text{f}(\text{t},\text{mn@boldsymbol}\theta )+\text{mn@boldsymbol}ϵ$$

$$b_{\hat{\theta }}(\text{mn@boldsymbol}{\theta }_t)=\mathbb{E}\{\hat{\theta }(\text{y})\}-\text{mn@boldsymbol}{\theta }_t$$

$$\mathrm{MSE}(\hat{\theta }):=\mathbb{E}\{(\hat{\theta }(\text{y})-\text{mn@boldsymbol}{\theta }_t{)}^2\}=\mathrm{Var}\{\hat{\theta }(\text{y})\}+b_{\hat{\theta }}(\text{mn@boldsymbol}{\theta }_t{)}^2.$$

$$f(t,e,K,P,\omega ,M_0)$$

$$\cos \nu$$

$$\sin \nu$$

$$E-e\sin E$$

$$p(\text{y}|\text{mn@boldsymbol}\theta ,\text{mn@boldsymbol}\beta )=\frac{1}{\sqrt{(2\pi {)}^N|\mathbf{V}(\text{mn@boldsymbol}\beta )|}}{e}^{-\frac{1}{2}(\text{y}-\text{f}(\text{t},\text{mn@boldsymbol}\theta ){)}^T\mathbf{V}(\text{mn@boldsymbol}\beta {)}^{-1}(\text{y}-\text{f}(\text{t},\text{mn@boldsymbol}\theta ))}$$

$${\hat{\theta }}_{\mathrm{ML}}=\arg \underset{\text{mn@boldsymbol}\theta \in \Theta ,\text{mn@boldsymbol}\beta \in B}{\max }p(\text{y}|\text{mn@boldsymbol}\theta ,\text{mn@boldsymbol}\beta ).$$

$${\hat{\theta }}_{LS}(\text{y})=\arg \underset{\text{mn@boldsymbol}\theta \in \Theta }{\min }(\text{y}-\text{f}(\text{mn@boldsymbol}\theta ){)}^T{\text{V}}^{-1}(\text{y}-\text{f}(\text{mn@boldsymbol}\theta )).$$

$$S:=\frac{1}{{\sigma }_k}=\frac{K_t}{\sigma }\sqrt{\frac{N-p}{2}}$$

$$b(e_t,S)\approx \frac{1}{S}\sqrt{\frac{\pi }{2}}{L}_{1/2}\left(-\frac{S^2{e}_t^2}{2}\right)-e_t$$

$$b(0,S)\approx \sqrt{\frac{\pi }{2}}\frac{1}{S}=\frac{\sigma }{K_t}\sqrt{\frac{\pi }{N-p}}=\sqrt{\frac{\pi }{4-\pi }}{\sigma }_{\hat{e}}.$$

$$L{R}_e:=\frac{\underset{\text{mn@boldsymbol}\theta \in {\Theta }_e}{\max }p(\text{y}|\text{mn@boldsymbol}\theta )}{\underset{\text{mn@boldsymbol}\theta \in \Theta }{\max }p(\text{y}|\text{mn@boldsymbol}\theta )}⩽e^{-\frac{1}{2}\beta }\text{,where }$$

$$\beta :=F_{{\chi }_{\rho }^2}^{-1}(1-\alpha )$$

$$\rho :=2+2S^{\prime 2}\frac{e^2}{1+e^2}-\frac{\pi e}{1+e^2}{L}_{\frac{1}{2}}\left(-\frac{S^{\prime 2}}{2}\right)L_{\frac{1}{2}}\left(-\frac{e^2{S}^{\prime 2}}{2}\right).$$

$$\mathrm{FAP}(e)=1-F_{{\chi }_{\rho }}^2(-2\ln L{R}_e).$$

$${\hat{\theta }}_{\mathrm{MAP}}=\arg \underset{\text{mn@boldsymbol}\theta \in \Theta }{\max }p(\text{mn@boldsymbol}\theta |\text{y}).$$

$$p(e|\text{y})$$

$$\begin{array}{rl}{e}_{\max }=\arg \underset{e\in [0,1]}{\max }p(e|\text{y})& ;e_{\mathrm{mean}}=\mathrm{mean}(p(e|\text{y}))\\ e_{\mathrm{med}}=& \mathrm{median}(p(e|\text{y}))\end{array}$$

$$\mathrm{MAE}:=\mathbb{E}\{|\hat{e}-e_t|\}.$$

$$R=\frac{\mathrm{Pr}\{e\ne 0|\text{y}\}}{\mathrm{Pr}\{e=0|\text{y}\}}=\frac{\mathrm{Pr}\{\text{y}|e\ne 0\}}{\mathrm{Pr}\{\text{y}|e=0\}}\frac{\mathrm{Pr}\{e\ne 0\}}{\mathrm{Pr}\{e=0\}}$$

$$\mathrm{Pr}\{e\in C|\text{y}\}={\int }_{C}p(e|\text{y})de=\alpha$$

$$p(\text{y}|\text{mn@boldsymbol}\theta ,{\sigma }_J)=\frac{1}{\sqrt{(2\pi {)}^N|{\mathbf{V}}_0+{\sigma }_J^2\mathbf{I}|}}{e}^{-\frac{1}{2}(\text{y}-\text{f}(\text{mn@boldsymbol}\theta ){)}^T({\mathbf{V}}_0+{\sigma }_J^2\mathbf{I}{)}^{-1}(\text{y}-\text{f}(\text{mn@boldsymbol}\theta ))}$$

$${\kappa }_{i,j}={\alpha }^2\exp \left[-\frac{1}{2}\left\{\frac{{\sin }^2[\pi (t_i-t_j)/\tau ]}{{\lambda }_p^2}+\frac{(t_i-t_j{)}^2}{{\lambda }_e^2}\right\}\right],$$

$$R_{\text{mn@boldsymbol}ϵ}(\text{M})$$

$$Q_{\text{mn@boldsymbol}ϵ}(\text{M})$$

$$b(0,S)\approx \sqrt{\frac{\pi }{2}}\frac{q}{S}.$$

$$a$$

$$x$$

$$p^{\prime }(e,\widetilde{\text{mn@boldsymbol}\theta })$$

$$z$$

$$I_E$$

$$\begin{array}{rl}\mathrm{Pr}\{e\in C|\text{y},p^{\prime }\}& ={\int }_C{\int }_{\widetilde{\text{mn@boldsymbol}\theta }}{p}^{\prime }(e|\text{y})\mathrm{d}\widetilde{\text{mn@boldsymbol}\theta }\mathrm{d}e\\ & =\frac{x{I}_{C\cap D}+z{I}_{C\cap \bar{D}}}{x{I}_D+z{I}_{\bar{D}}}\end{array}$$

$$\widehat{\mathrm{Pr}}\{e\in C|\text{y},p^{\prime }\}=\frac{x|C\cap D|+z|C\cap \bar{D}|}{x|D|+z|\bar{D}|}.$$

Full text

††pubyear: 2019††pagerange: Bias and robustness of eccentricity estimates from radial velocity data–LABEL:lastpage

1 Introduction

The nearly coplanar and circular orbits of the Solar system have long been an argument in favour of Laplace and Kant’s theory of formation of planets in a disk (Swedenborg, 1734; Kant, 1755; Laplace, 1796). The first observations of exoplanets suggested that such low eccentricities were rather the exception than the rule. The “eccentricity problem”, along with the need to envision migration scenarios for hot Jupiters, triggered several theoretical studies which explored migration scenarios after the disk has dissipated. The prediction of these models were compared to measured eccentricities. For instance Jurić & Tremaine (2008) and Petrovich et al. (2014) evaluate the likelihood of formation scenarios of hot and warm Jupiters through their agreement with observed eccentricity distributions.

For a radial velocity data set on a given star, one wants to extract two pieces of information about the eccentricity. First, a best candidate value (the estimation problem) and what are the eccentricity values that are incompatible with the data (the hypothesis testing problem). It is in particular interesting to test if an eccentricity is non zero. Both problems do not have completely obvious solutions. For instance it is known since Lucy & Sweeney (1971) that when the true eccentricity is small, its least square estimate is biased upwards. Other aspects of the estimation and hypothesis testing problems have been addressed in the exoplanet community.

In Shen & Turner (2008) the bias was found to depend on the signal to noise ratio (denoted by SNR) as well as on the time span of the observations. This was confirmed by Zakamska et al. (2011), which further showed that the bias depends on the phase coverage, and updated the Lucy & Sweeney (1971) null hypothesis test to determine if a null eccentricity can be rejected or not. They also propose metrics for evaluating the quality of a data set. O’Toole et al. (2009a) showed that error bars on eccentricity from least square are under-estimated by a factor 5-10. Brown (2017) shows that there might be orbits with very high eccentricities with similar goodness-of-fit as a low eccentricity one. Pont et al. (2011), Husnoo et al. (2011) and Husnoo et al. (2012) used Bayesian Information Criterion to confirm non-zero eccentricities. More recently, Bonomo et al. (2017b, a) assessed the evidence in favour of eccentric solutions with Bayes factors. A Bayesian test with a physically motivated prior on eccentricity was devised by Lucy (2013). Also, Anglada-Escudé et al. (2010), Wittenmyer et al. (2013) and Kürster et al. (2015) note that two planets in 2:1 mean motion resonances can appear as an eccentric planet, and propose ways to disentangle those cases. This problem has also been studied in Boisvert et al. (2018); Nagel et al. (2019); Wittenmyer et al. (2019).

The fact that eccentricity estimates can be spuriously high for a given planet gives reasons for concern on the eccentricity distributions. The estimation of those has been tackled in Hogg et al. (2010), which computes the posterior of the eccentricity probability distribution itself.  Zakamska et al. (2011) consider the accuracy of the eccentricity catalogues obtained by Bayesian point estimates. They show that for single planet populations contaminated by white noise, estimating the eccentricity via the maximum of the marginalised posterior distribution of eccentricity with a free jitter term gives satisfactory retrieval of the input population. Furthermore, it has been noted by Cumming (2004) that high eccentricity orbits $\gtrsim 0.6$ are more difficult to detect at fixed semi-amplitude. For a fixed mass, the detection bias is less strong (Shen & Turner, 2008).

We contribute to this series of work by studying in depth the bias and robustness of eccentricity estimates from radial velocity data (note that our analysis also applies to astrometric measurements and to estimates of semi-amplitude and inclination). We proceed in two steps. First, we highlight key properties of classical estimators, in order to have a consistent view of eccentricity estimation. The following questions are then considered: is the eccentricity inference robust to modelling errors? By that, we refer to wrong noise models, planetary companion too small to be detected, etc. If not, how to mitigate the problem? One could encounter a situation similar to the spectroscopic binaries, where proximity effects or gas streams lead to spurious high-significance eccentricities if not properly accounted for (Lucy, 2005).

The article is structured as follows. In section 2 we study the behaviour of eccentricity estimates when the model is correct. The least square estimates as well as Bayesian ones are studied, and is is shown that the latter are less biased at low eccentricity. The problem of spurious local $\chi^{2}$ minima at high eccentricity is also tackled, in particular through the Proxima b case (Anglada-Escudé et al., 2016). Section 3 is devoted to studying the robustness of the estimates when the numerical method, the model or the prior is poorly chosen. Finally in section 4, we consider ways to check the validity of a noise model. Our methods are illustrated with the 55 Cnc HIRES data set in section 5. In section 6, we conclude with a step-by-step procedure to obtain reliable eccentricities and present perspectives for future work .

2 Point and interval eccentricity estimates

2.1 Problem definition

2.1.1 Point and interval estimates

Let us first define the problem under study precisely. Some generic symbols, used throughout the text, are summarized in table 1.

Let us consider a time series of $N$ observations, modelled as a vector $\textbf{{y}}=(y(t_{k}))_{k=1..N}$, such that

[TABLE]

where $\textbf{{t}}=(t_{k})_{k=1..N}$ is the vector of measurement times, $f$ is a deterministic model depending on parameters $\mn@boldsymbol{\theta}\in\mathbb{R}^{p}$ and $\mn@boldsymbol{\epsilon}$ is a random variable modelling the noise. An estimator of $\mn@boldsymbol{\theta}$ is a function $\widehat{\theta}$ of the data $\textbf{{y}}(\textbf{{t}})$ whose output is wanted to be close to the true value of $\mn@boldsymbol{\theta}$, denoted by $\mn@boldsymbol{\theta}_{t}$, in a sense chosen by the data analyst. If the mean value of $\widehat{\theta}(\textbf{{y}})$ ($\mathbb{E}\{\widehat{\theta}(\textbf{{y}})\}$) is not equal to $\mn@boldsymbol{\theta}_{t}$, the estimator is said to be biased and

[TABLE]

is called the bias of the estimator $\widehat{\theta}$ in $\mn@boldsymbol{\theta}_{t}$. A common metric for the accuracy of an estimator is the mean squared error (MSE), linked to bias and the variance of the estimator, $\mathrm{Var}\{\widehat{\theta}(\textbf{{y}})\}$, via

[TABLE]

The other problem we consider is to have a testing procedure to reject or not certain values of the eccentricity. We are now interested in rejecting the hypothesis that $e\in C$ where $C$ is a subset of $[0,1]$. More precise definitions are given in the relevant sections (2.2.4 and 2.3.2).

In the present section 2, we describe tools for the estimation and hypothesis testing problems and present some of their properties. In section 3, we study the reliability of these tools when the model is incorrect. By convention, in the following, radial velocity signals are in m.s-1. The analysis is unchanged for other units as long as the signal to noise ratio is identical.

2.1.2 Model

The concern of the present work is the estimation of eccentricity from radial velocity data. The model of a radial velocity planetary signal is recalled below

[TABLE]

The symbols $t,e,K,P,\omega,M_{0}$ designate respectively the measurement time, eccentricity, semi-amplitude, period, argument of periastron and mean anomaly at $t=0$. The symbols $E$ and $\nu$ denote the eccentric and true anomalies.

We assume a Gaussian noise model, such that the likelihood function is

[TABLE]

where $\textbf{{f}}(\textbf{{t}},\mn@boldsymbol{\theta})$ is a sum of Keplerian functions defined as (4) possibly plus some other model features (offset, trend…). The covariance matrix $\mathbf{V}$ is parametrized by $\mn@boldsymbol{\beta}$ and the suffix $T$ denotes the matrix transposition. The explicit expression of $\textbf{{f}}(\textbf{{t}},\mn@boldsymbol{\theta})$ and $\mathbf{V}(\mn@boldsymbol{\beta})$ will be given in the relevant sections.

The features of least square and Bayesian estimates are now studied respectively in sections 2.2 and 2.3.

2.2 Least square estimate

2.2.1 Bias of the non linear least square

A common parameter estimator is the maximum likelihood $\widehat{\theta}_{\mathrm{ML}}$. For the model of eq. (8),

[TABLE]

When the parameters of the covariance, $\mn@boldsymbol{\beta}$, are fixed, maximising the likelihood comes down to the least square problem,

[TABLE]

When the model $\textbf{{f}}(\mn@boldsymbol{\theta})$ depends linearly on $\mn@boldsymbol{\theta}$, the least square estimate is unbiased. This is in general untrue when $\textbf{{f}}(\mn@boldsymbol{\theta})$ is non linear, which has been noted for instance by Hurwicz (1950) and discussed in Hartley (1964); Bates & Watts (1980); Cook & Witmer (1985); Firth (1993). The model we are concerned with (eq. (4) - (7)) is non linear, and indeed eccentricities obtained by least square are biased.

2.2.2 Bias and uncertainty at low eccentricity

We begin with a numerical experiment. We generate Keplerian signals of eccentricity 0.06, 0.5 and 0.9 with fixed $\omega$, $M_{0}$ and $K=3$ m.s-1. The array of time t is the 30 first measurements of GJ 876 (Correia et al., 2010). We generate 1000 realisations of a white Gaussian noise with standard deviation of 1 m.s-1. For each realisation of the noise, a non-linear Keplerian model and a constant are fitted with a Levenberg-Marquardt algorithm. The minimisation step is scaled so that the eccentricity never exceeds 0.999. The values of the estimates ($\widehat{k}:=\widehat{e\cos\omega},\widehat{h}:=\widehat{e\sin\omega}$) are reported as crosses (red, purple and blue resp. for $e=0.06$, 0.5 and 0.9) in Fig. 1. The distributions of $\widehat{k}$ and $\widehat{h}$ are fairly isotropic for $e_{t}=0.06$ and $e_{t}=0.5$. For $e=0.9$ there seems to be more complicated phenomena at work. The distribution has no circular symmetry and in some cases $\widehat{e}$ is stuck at its maximal value, 0.999 (see section 2.2.3).

In appendix A, it is shown that provided $e$ is small enough ($\leqslant 0.2$) and the number of observation is sufficient, $\widehat{k}$ and $\widehat{h}$, follow independent Gaussian laws of same variance. This property allows us to understand the bias qualitatively. In Fig. 2, the distribution of eccentricity estimates is represented. The pair of true values $k_{t},h_{t}$ is represented by a yellow star. The bold black line delimits the region where 99$\%$ of the estimates are located. When the estimate falls in the blue-coloured region, the eccentricity is under-estimated. When it falls in the red-coloured region, the eccentricity is over-estimated. As the volume of higher eccentricity models is larger in the vicinity of $k_{t},h_{t}$, the eccentricity is more probably over-estimated. Informally, there are more and more models with eccentricity $e$ as $e$ grows.

Furthermore, it is possible to obtain an analytical approximation of the bias. Since $\widehat{k}$ and $\widehat{h}$ approximately follow a joint Gaussian distribution with same variances, $\widehat{e}=(\widehat{k}^{2}+\widehat{h}^{2})^{1/2}$ follows a Rice distribution, as noted in Shen & Turner (2008). Interestingly enough, the Rice distribution appears as a very good model for the eccentricity densities of the inner planets of the Solar System, resulting from chaotic diffusion (Laskar, 2008). Supposing that the measurement noise is white with standard deviation $\sigma$ and that $p$ parameters are fitted, within our approximation, $\widehat{k}$ and $\widehat{h}$ have a standard deviation $\sigma_{k}=\sigma_{h}=\sigma/K_{t}\sqrt{2/(N-p)}$. Defining the SNR $S$ as

[TABLE]

and denoting by $e_{t}$ the true eccentricity, the bias is

[TABLE]

where $L_{1/2}$ is the Laguerre polynomial of order $1/2$. When $e_{t}=0$, the eccentricity follows a Rayleigh distribution and eq. (12) reduces to a very simple expression,

[TABLE]

Eq. (13) is identical to equation (18) of Lucy & Sweeney (1971) except that we are able to derive the effect of the correlations between parameters on the SNR through the term $-p$ (see Appendix A for justification).

Formula (13) is useful to see a few trends: the bias is proportional to the uncertainty on $k$ and $h$, which is proportional to the inverse of $K_{t}$ and $\sqrt{N-p}$. As a consequence, the bias increases as the SNR decreases, i.e. as $\sigma$ increases or as $K_{t}$ or $N$ decrease. This is also found by simulations in Shen & Turner (2008) and Zakamska et al. (2011) for Bayesian estimates. We add that increasing the number of fitted parameters $p$, increases the bias.

There are particular cases where the correlations between parameters drastically increase the uncertainties on $k$ and $h$ and therefore increase the bias, so that formula (12) should be taken as a lower bound.

However, the fact that the bias is approximately proportional to the standard deviation with a factor $\sqrt{\pi/(4-\pi)}$, as in eq. (13) stays true. This fact is remarkable because it means that the accuracy of the estimate (seen as the mean squared error (3)) is proportional to its precision (seen as the standard deviation $\sigma_{\widehat{e}}$). In appendix A.2, we show that poor phase coverage or short observational baseline affect the accuracy of the eccentricity estimate insofar as they decrease the precision of the estimate.

2.2.3 Local minima at high eccentricities

As shown in Baluev (2015), the number of local $\chi^{2}$ minima increases significantly in the high eccentricity region. These minima might lead a local minimisation algorithm or a Monte Carlo Markov Chain (MCMC) to be stuck in the wrong region of the parameter space. We here aim at quantifying and understanding this feature. In this section, we provide a summary of our results, whose precise study is in appendix B. These results are:

• •

An incorrect estimation of the noise level can lead to spurious deep local minima at high eccentricities.

• •

As the SNR degrades, the probability of missing the global minimum by a non linear fit initialized on a circular orbit increases.

• •

There is a geometrical interpretation of the numerous local minima at high eccentricity: the set of models with fixed eccentricity explore more and more dimensions of the $N$-dimensional sample space as eccentricity grows.

Let us illustrate the first point on Proxima b. Brown (2017) re-analyses the data of Proxima b, a $\approx 1.27$ $\mathrm{M}_{\oplus}$ planet orbiting the M star Proxima Centauri with a period of 11.186 days (Anglada-Escudé et al., 2016), and finds that there are local minima at eccentricity 0.75 and 0.95, the 0.95 eccentricity being the global least square fit.

We compute a Keplerian periodogram (O’Toole et al., 2009b; Zechmeister & Kürster, 2009; Baluev, 2015) in the vicinity of the 11.186 period. That is Keplerian models (eq. (4) - (7)) are fitted for a grid of periods, argument of periastron and eccentricity. We then represent per eccentricity the $\chi^{2}$ minimised over all other parameters. Using the nominal uncertainties, a single Keplerian model plus one offset per instrument, a linear and quadratic trend, we obtain the red curve in Fig. 3.

The curve displays three local minima, the deepest being at eccentricity 0.92. However, let us note that the best fit gives $\chi^{2}$ of 1057, while there are only 214 measurements. This gives a reduced $\chi^{2}$ of 5.16, which is unrealistic. We here simply add a constant jitter term in quadrature with the nominal error bars to obtain a reduced $\chi^{2}$ of 1 at the best fit. The minimum $\chi^{2}$ as a function of eccentricity so obtained is represented in blue in Fig. 3. The global minimum now occurs at $e=0.17$. We interpret the global minima at eccentricity 0.92 as an artefact of an incorrect estimation of the error bars.

The same Keplerian periodogram calculations can be done on simulated data sets with different noise levels and different numbers of measurements. We simulate such systems with eccentricity sampled from a uniform distribution on $[0,0.999]$, and count how many of them that have a SNR between 0 and 5, 5 and 10 etc. have $k=1,2,3...$ local minima. The histogram of Fig. 4 is obtained. Also, for each bin of SNR, we compute the proportion of systems where the global minimum is not the closest to 0 (as in the case on Fig. 3, red curve) and therefore a local minimisation should miss the global minimum. It appears that as the SNR increases, the fraction of cases where the global minimum is missed is decreasing, though not reaching zero.

Let us now briefly comment on the geometrical interpretation of the higher number of local minima at high eccentricities. Finding the best fitting model amounts to finding the closest model to the observation in a geometrical sense. We consider the figure drawn in $\mathbb{R}^{N}$ by all the models that have an eccentricity $e$ and a period $P$, denoted by $\mathcal{M}_{e,P}$. This figure might explore more or less dimensions. For instance, if it is close to a plane, it is nearly confined to a two-dimensional space. Otherwise, exploring many dimensions traduces a “rough” surface, which increases the chances of finding a local minimum of distance to the data. By a procedure based on singular value decomposition, detailed in appendix B and in Hara 2017 (PhD thesis), it is possible to obtain an approximate number of dimensions explored by $\mathcal{M}_{e,P}$ as a function of $e$. Table 2 shows such calculations for the 214 measurement times of GJ 876 (Correia et al., 2010). As eccentricity increases, $\mathcal{M}_{e,P}$ explores more and more dimensions.

2.2.4 Interval estimates

As said in section 1, one does not only want to obtain a value of the eccentricity with error bars, but also to test whether a certain value of the eccentricity is compatible with the data. This can be done in the frequentist setting with interval estimates. Since the focus will be put on Bayesian estimates in section 3, we simply here give their definition and refer the reader to appendix C for their derivation and detailed study.

The hypothesis that the eccentricity has a certain value $e$ is rejected with a confidence level $\alpha$ if the likelihood ratio $LR_{e}$ satisfies

[TABLE]

where $\Theta_{e}$ is the set of parameters that have all eccentricity $e$, $p(\textbf{{y}}|\mn@boldsymbol{\theta})$ is the likelihood, $F_{\chi^{2}_{\rho}}^{-1}$ is the inverse cumulative distribution function of a $\chi^{2}$ law with $\rho$ degrees of freedom, $S^{\prime}=(\sigma/K_{t})\sqrt{2/N}$ and $L_{\frac{1}{2}}$ is the Laguerre polynomial of order $1/2$. We also fit a free jitter term so that the reduced $\chi^{2}$ equals one. Conversely, for a certain measured $LR_{e}$ the FAP$(e)$ of the hypothesis $e_{t}=e$ is defined as

[TABLE]

The confidence interval of eccentricity is the set of $e$ with $\mathrm{FAP}(e)$ greater than a certain threshold (for instance $0.05$).

2.3 Posterior distributions

2.3.1 Point estimators

The previous sections are devoted to the study the least square eccentricity estimate. However, the standard practice in the exoplanet community is rather to compute the posterior probability $p(\mn@boldsymbol{\theta}|\textbf{{y}})=p(\textbf{{y}}|\mn@boldsymbol{\theta})p(\mn@boldsymbol{\theta})/p(\textbf{{y}})$ of the orbital elements $\mn@boldsymbol{\theta}$ using Monte Carlo Markov Chains algorithms (e.g. Ford, 2005, 2006).

From such posterior distributions, one can compute the orbital elements corresponding to the maximum a posteriori (MAP)

[TABLE]

For a simple Keplerian fit $\widehat{\theta}_{\mathrm{MAP}}=(K,e,P,\omega,M_{0})^{\mathrm{MAP}}$, an eccentricity estimate is then obtained by $e^{\mathrm{MAP}}$. Alternately, one can compute the marginal distribution

[TABLE]

and its maximum, mean and median

[TABLE]

where $\tilde{\mn@boldsymbol{\theta}}$ is the vector of parameters that are not eccentricity, and $\Theta_{e}$ is the space of parameters where the eccentricity is equal to $e$. A standard result is that $e_{\mathrm{mean}}$ is the estimator that minimizes the mean squared error (see eq. (3)). Also, $e_{\mathrm{med}}$ minimizes the mean absolute error (MAE), defined as

[TABLE]

The estimators (20) are in general less than the maximum likelihood or the maximum a posteriori. This is shown with a numerical experiment. We generate a circular planet of semi-amplitude 3.5 m/s and 100 realizations of Gaussian white noise at 2 m/s. The estimates $e_{\mathrm{ML}}$ and $e_{\mathrm{mean}}$ (eq. (20)) are computed with a Monte-Carlo Markov Chain (MCMC) algorithm. The sampler is the same as in Delisle et al. (2018), based on the adaptive Metropolis sampler of Haario et al. (2001). The model consists of a Keplerian, an offset and a free jitter term, with uniform priors on all parameters. Fig. 5 shows experimental distribution function of $e_{\mathrm{ML}}-e_{\mathrm{mean}}$. The condition $e_{\mathrm{ML}}>e_{\mathrm{mean}}$ is verified in 79 cases out of 100, with a mean value of $e_{\mathrm{ML}}-e_{\mathrm{mean}}$ equal to 0.0944.

The efficiency of the estimates (20) is understandable in terms of trade-off between model simplicity and agreement with the data. In section 2.2.2 and 2.2.3, it appeared that for eccentricities $e_{t}$<0.2, in the vicinity of the true model, there is a larger volume of models with eccentricities $e>e_{t}$ than $e<e_{t}$. The integration over the domain $\Theta_{e}$ in eq. (19) penalizes the models with high complexity, which here are the high eccentricity models. This is comparable to the penalization of models with too many planets by the marginal likelihood in the context of exoplanets detection (e.g. Nelson et al., 2018). In section 2.3.2, we argue in favour of reporting $e_{\mathrm{mean}}$ and/or $e_{\mathrm{med}}$, as $e_{\mathrm{max}}$ is too biased towards low eccentricities.

2.3.2 Hypothesis testing

Several of the works cited in section 1 address the question of whether an eccentricity should be set to zero or not, which is a model selection problem. It can be addressed by computing the ratio of posterior likelihood, or odds ratio, of the models (e.g. Kass & Raftery, 1995),

[TABLE]

where $\mathrm{Pr}\{\textbf{{y}}|e\neq 0\}=\int_{\mn@boldsymbol{\theta}\in\Theta}p(\textbf{{y}}|\mn@boldsymbol{\theta})p(\mn@boldsymbol{\theta})d\mn@boldsymbol{\theta}$. If this ratio is greater than a certain value, then one favours $e\neq 0$ over $e=0$. This methodology has been used in Bonomo et al. (2017b, a). As the number of samples $N$ tends to infinity, assuming $\mathrm{Pr}\{e\neq 0\}=\mathrm{Pr}\{e=0\}$, the odds ratio is equivalent to the Bayesian Information Criterion (BIC, Schwarz, 1978), as used by Pont et al. (2011); Husnoo et al. (2011, 2012).

More generally, one can compute a credible set, that is a set of $e$, denoted by $C\subset[0,1]$ such that

[TABLE]

where $\alpha\in[0,1]$ is a probability. The set $C$ is in general taken as an interval but this need not be the case. Let us note that this approach is also a ratio of posterior likelihood as in Eq (22). If $e=0$ is given a non null probability, then the prior probability takes the form $p(e)=p(0)\delta(0)+p(e)$ where $\delta$ is the Dirac function. We now focus on credible intervals, (eq. (23)) since these are more widely used and give finer information.

3 Robustness of eccentricity estimates

3.1 Problem statement

In section 2, several tools to make inferences on eccentricity were presented. We now study whether these are reliable even if the adopted model is incorrect. What we call a model is a couple of prior and likelihood functions. We assume that the orbital elements are distributed according to a true prior $p_{t}(\mn@boldsymbol{\theta})$ and the observations have a true likelihood $p_{t}(\textbf{{y}}|\mn@boldsymbol{\theta})$ ($\mathcal{M}_{t}$) while the analysis is made with the model $\mathcal{M}$: $p(\mn@boldsymbol{\theta}),p(\textbf{{y}}|\mn@boldsymbol{\theta})$. The model $\mathcal{M}$ can be too simple: missed planetary signal, non modelled correlated noise or too complicated: for instance Gaussian processes are known to be very flexible, possibly too much.

The case where the model is too complicated will not be treated in detail. We simply point out that from section 2.2, we expect that the bias is higher than for a correct model. Indeed, as the model grows in complexity the correlation between parameters increases, therefore the error on eccentricity $\sigma_{e}$ increases. This might be problematic as the bias is proportional to $\sigma_{e}$ (see eq. (13). However, the error bars broaden, so that having too complex a model is unlikely to produce spurious conclusions. On the contrary, as we shall see, simplistic models can be problematic.

In the following sections, we study the effect of the noise level estimate (3.2), numerical effects (3.3), incorrect noise models (3.4), priors (3.5) and the comparison of two models: one eccentric planet or two planets in 2:1 mean motion resonance (3.6).

3.2 Noise level

In section 2.2.2 and 2.2.3, it appeared that an incorrect estimate of the noise norm leads to an underestimated bias and to spurious local minima at high eccentricity. As a consequence, it is key, as is standard practice, to adjust at least an extra jitter term $\sigma_{J}$ in the likelihood,

[TABLE]

where $\mathbf{V_{0}}$ is the nominal covariance, $\mathbf{I}$ the identity matrix and $\textbf{{f}}(\mn@boldsymbol{\theta})$ is the signal model, containing Keplerians and possibly other features. All the following analyses are made with the model (24).

3.3 Numerical effects

It appeared in section 2.2.2 that uncertainties on the eccentricity estimates increase the biases. This is also valid for the uncertainties stemming from the numerical methods used to compute the orbital elements.

As noted by Eastman et al. (2013), there is a specific error in the implementation of the Metropolis-Hastings algorithm that worsens the bias, when the true eccentricity is close to zero and when the parameter space is parametrized by $(e,\omega)$ instead of $(k,h)$ or $(\tilde{k},\tilde{h})$. This error consists in not recording the value of a proposed parameter in the chain if it is rejected.

More generally, credible intervals or Bayes factors can be unreliable if the numerical schemes have not converged.  Hogg & Foreman-Mackey (2018) gives several ways to check for convergence of MCMC algorithms. In the following analyses, the convergence diagnostic is the number of effective samples, as computed in Appendix A of Delisle et al. (2018) (see also Sokal (1997)). This number, $N_{\mathrm{eff}}$ is interpretable as an equivalent number of independent samples from the posterior distribution. Then, for instance, the numerical uncertainty on the mean of the marginal posterior $p(\theta_{0}|\textbf{{y}})$ of a parameter $\theta_{0}$ scales as $\sigma_{p}/\sqrt{N_{\mathrm{eff}}}$, where $\sigma_{p}$ is the standard deviation of $p(\theta_{0}|\textbf{{y}})$.

The first result we show is that at low SNR, the convergence is slower. This is likely due to the existence of local minima at high eccentricity (see fig. 4, section 2.2.3). We simulate a one planet system on CoRoT-9 28 measurements (Bonomo et al., 2017b), the eccentricity is generated with a Beta distribution (a=0.867, b=3.03), angles are uniformly distributed and the period is fixed to 95 days. The nominal measurement errors are normalised to obtain a mean variance of 1. A Gaussian white noise following the normalised nominal errors plus a one meter per second white noise is added. A hundred data sets with orbital elements and noise sampled according to their distributions are created. This is done for $K=2.5$ m.s-1 and $K=5$ m.s-1, which corresponds to SNR 4.2 and 8.5 (as defined in eq. (11)), so to obtain $2\times 100=200$ data sets. Finally, a MCMC is performed on each of them with model (24) with 1,100,000 samples and a burn-in phase of a fourth of the total number of samples. The average number of independent samples is 7300 vs 12000, and the chain do no not reach a 1000 efficient samples in 28 and 5 cases respectively. In conclusion, low SNR signals should be treated with particular attention since the number of independent samples of the MCMC is smaller, while chains are initialized with a high eccentricity due to the higher bias of least square estimates at low SNR (see eq. (13)).

As an example of the problem associated with not having enough samples, this experiment is repeated with an injected eccentricity of 0. We consider that the “small eccentricity” hypothesis is rejected if $p(e\in[0,0.05])<0.05$. Let us consider the samples obtained for the experiment described above with $K=2.5$ m.s-1. Depending on whether we take the first fourth of the samples or all of them, the fraction of rejection of the small eccentricity hypothesis goes respectively from 12 to 4 $\%$.

In the following sections, simulations are taken into consideration if they reach an effective number of independent samples greater than 1000, so that the mean of the eccentricity posterior (of variance $\sigma_{e}$) is known with at least a $\approx\sigma_{e}/\sqrt{N_{\mathrm{eff}}}\approx\sigma_{e}\times 3\%$ accuracy. Since we compare different types of noise, we do not require an extremely good precision on the posterior distribution (also the following result stay approximately identical if we take runs with at least 2000, 5000 efficient samples). In practice, it is safer to ensure that 10,000 independent samples are reached to obtain a $\sigma_{e}/\sqrt{10000}=\sigma_{e}\times 1\%$ accuracy on the posterior mean.

3.4 Incorrect noise model

3.4.1 Non Gaussian noise

The first question we address is whether, when the noise is non Gaussian, using the model (24) leads to spurious inferences on eccentricity. For the sake of brevity, we here only report the results of the analysis done in appendix D. The non Gaussianity of the noise can only lead to slightly underestimating or overestimating the error bars, when the noise distribution is very heavy-tailed. We found that the estimates of the eccentricity and the error bars are mostly sensitive to the covariance of the noise and the following sections are focused on this aspect.

3.4.2 Incorrect covariance: simulations

To study the effect of the noise covariance on eccentricity estimates, we proceed as follows. We consider the 28 measurement times of CoRoT-9, spanning on 1700 days, and generate a circular planet at 95 days, denoted by $y_{\mathrm{planet}}(t)$. We then generate 100 realisations of white noise, $y_{\mathrm{noise}}(t)$, and for each of them, a non modelled noise is added. Six different types of such noises, described below, are considered. In total, 6$\times$100 = 600 data sets are obtained. The signal generated is then of the form $y(t)=y_{\mathrm{planet}}(t)+y_{\mathrm{noise}}(t)+y_{k}(t)$ with $k=0..5$. $y_{\mathrm{noise}}(t)$ is generated according to the nominal uncertainties, which are CoRoT-9 uncertainties scaled so that their mean square is equal to 1 m.s-1. Such signals are generated in three contexts: with a semi-amplitude of the planet at 95 days of $2.5$, $3.5$ or $5$ m.s-1. The methodology described aims at evaluating if, when the noise model is incorrect, circular planets tend to appear as eccentric. This experiment is also done with the eccentricity of the planet drawn from the same Beta distribution as the prior, to evaluate the impact of the noise on the estimates. The simulation where the eccentricity is fixed to zero and the one where it follows a beta distribution are respectively called $S_{0}$ and $S_{e}$.

In each simulation, on each of the 3$\times$600 data sets, the posterior distribution of the orbital elements is computed using the model (24) and priors given in table 5. $1,100,000$ samples are computed, the first fourth being the burn-in phase. The algorithm is an adaptive Metropolis algorithm as in Delisle et al. (2018), and the convergence is checked by calculating the effective number of independent samples (see Delisle et al. (2018), Appendix A).

Since our goal is to test the effect of the noise nature, and not its level, we impose that the norm of the non modelled noise is such that $\|\textbf{{y}}_{0}\|=\|\textbf{{y}}_{1}\|=...\|\textbf{{y}}_{5}\|=\gamma$. For each of the 100 realisations of white noises, we draw $\gamma$ from a $\chi^{2}$ law with $N=28$ degrees of freedom. The $y_{k}$ are defined as follows:

• $y_{0}$

: white, Gaussian noise identically distributed.

• $y_{1}$

: also white Gaussian noise but with different variances. The variances are drawn from a Laplace law so to obtain a wide range of value.

• $y_{2}$

: A circular planet that is too small to be fitted, the period of the added planet $P$ is drawn from a log-normal distribution until a period is found such that $P$ differs from 95 days of at least 20%. This value is chosen to avoid the lowest probability region of period ratios of planet pairs found by Kepler (Steffen & Hwang, 2015).

• $y_{3}$

: planet in resonance with the injected planet, in 1:2 or 3:2 resonance, inner or outer with probability 1/2.

• $y_{4}$

: a Gaussian correlated noise with covariance $\kappa$

[TABLE]

as in Haywood et al. (2014). We use the values of the Evidence Challenge (Nelson et al., 2018), $\alpha=\sqrt{3}$ m/s, $\lambda_{e}=50.0$ days, $\lambda_{p}=0.5$ (unitless), and $\tau=20.0$ days.

• $y_{5}$

: same as $y_{4}$ but with values $\alpha=\sqrt{3}$ m/s, $\lambda_{e}=50.0$ days, $\lambda_{p}=0.3$ (unitless), and $\tau=30.0$ days.

The rationale behind taking $y_{1}$ as such is to emulate the effect of mild outliers, not obvious enough to be completely discarded. As the Laplace distribution has heavy tails, it generates values of the variances that are very different from each other. $y_{2}$ and $y_{3}$ are chosen as such because the strongest resonances found in Kepler data are the 3:2 and 1:2 (or close to) resonances. In first approximation, the real period ratio distribution is a combination of these two artificial distributions. Finally, $y_{4}$ and $y_{5}$ are two types of stellar noises, the second one having a slightly stronger periodic component in the covariance.

For the simulation where the eccentricity is fixed to zero, $S_{0}$, on each simulation we compute the probability $\mathrm{Pr}\{e\in[0,0.05]\}$. We report the number of simulations where this quantity is below a threshold $\alpha=0.05$ and that have at least a thousand effective samples. The results are shown on Fig. 6 for different values of the semi-amplitude of the input planet ($K=2.5$, 3.5 or 5 m.s-1, respectively blue, red and yellow markers). These values of $K$ were chosen to be close to relatively low SNRs. The small eccentricity rejection rate is approximately constant for the different noises except for a resonant companion, where the non zero eccentricity is, on average, rejected in 26% of the cases, that is 13 times more often than for the other noise models. A non resonant companion also might induce an increased rate of rejected zero eccentricity (4.2% on average vs 1.3% for white noise). The power of the stellar noise is stronger around 20 and 30 days, which is not on a harmonic of the period of the planet (95 days), and therefore they lead to an even lower rejection of the low eccentricity scenario (see section 3.4.3).

In the simulation $S_{e}$, we compute the absolute value of the difference between the estimated eccentricity and the true one for the three estimators (20). In all cases, we find that the median and mean of the posterior distribution of eccentricity are very similar and largely outperform the other estimates in terms of mean squared error (MSE) (eq. (3)) and mean absolute error (MAE) eq. (21)). We have seen in section 2.3 that for a correct model, the mean and median have respectively minimal MSE and MAE, this is indeed the case on the simulations with the white noise. The MSE and MAE over the 6 types of noise are reported in table 3. The estimator $e_{\mathrm{mean}}$ and $e_{\mathrm{med}}$ are more accurate when the noise model is incorrect both for MSE and MAE, as a consequence we deem them as the best ones overall. An argument in favour of $e_{\text{mean}}$ is that it has minimal MSE (eq. (3)). The MSE penalizes the square of the difference between the estimated and true value. Multiplying by two this difference multiplies by four the cost of the error. Therefore, estimators with small MSE are less likely to produce large errors.

This result seems in contradiction with Zakamska et al. (2011), which finds the mode of the posterior to be less biased. However, they consider cases where the eccentricity is small. We also find that for small eccentricities, $e_{\mathrm{max}}$ is less than $e_{\mathrm{med}}$ and $e_{\mathrm{mean}}$. Due to the Beta prior in eccentricity, it happens that the posterior is bi-modal and $e_{\mathrm{max}}=0$. When $e$ follows the prior distribution, $e_{\mathrm{med}}$ and $e_{\mathrm{mean}}$ are more accurate, at least in terms of MSE and MAE.

We now only present in Fig. 7 the performance in terms of MAE of the estimator $e_{\mathrm{mean}}$ (circles) and the maximum likelihood estimator $e_{\mathrm{ML}}$ (crosses), for comparison (MSE behaves similarly).

In all cases $e_{\mathrm{mean}}$ shows better performance. The MAE does not heavily depend on the type of noise. However, we do observe a slight increase of the error for the unseen resonant companion. Let us also note that as the signal semi-amplitude increases, the difference of performance between the maximum likelihood and $e_{\mathrm{mean}}$ becomes less clear.

3.4.3 Interpretation: noise power at the planet semi-period

In sections 3.4.1 we showed that the error on eccentricity is mainly determined by the true covariance of the noise. In section 3.4.2 however, the simulated stellar noises did not yield particularly high errors on the eccentricity. We now show that the property of correlated noises most impacting eccentricity estimates is their power at the semi-period of the planet of interest.

Let us consider a signal $\textbf{{y}}=\textbf{{y}}_{0}+\mn@boldsymbol{\epsilon}$, where $\textbf{{y}}_{0}$ is a circular orbit of period $P$, and $\mn@boldsymbol{\epsilon}$ is an unknown stochastic signal, that the data analyst supposes to have covariance V. Denoting by $\omega_{0}:=2\pi/P$, we define the $N\times 2$ matrix $\textbf{{M}}_{k}=[\cos(k\omega_{0}\textbf{{t}}),\sin(k\omega_{0}\textbf{{t}})]$ and $\textbf{{P}}:=\textbf{{M}}_{k}(\textbf{{M}}_{k}^{T}\textbf{{V}}^{-1}\textbf{{M}}_{k})^{-1}\textbf{{M}}_{k}^{T}\textbf{{V}}^{-1}$ the projection matrix onto the space spanned by the columns of $\textbf{{M}}_{k}$. More generally, for a projection matrix P onto a vector space M we define

[TABLE]

The rationale of defining these quantities is to identify if, assuming a covariance V a noise is more correlated to a certain space than to its orthogonal. If $\mn@boldsymbol{\epsilon}$ is a white noise, then it is not particularly correlated to any particular space, so that $R=1$. In the limit cases where $\mn@boldsymbol{\epsilon}$ lies in, or is orthogonal to the space spanned by $\textbf{{M}}_{k}$, then $R=N/2$ resp. 0.

To test the influence of $Q_{\mn@boldsymbol{\epsilon}}(\textbf{{M}}_{2})$ on the eccentricity error we proceed as follows. We consider an array of measurement times from a real system and generate a circular signal $y_{0}$ plus a noise with a certain true covariance. In order to obtain noises with very different spectral contents, we proceed as follows. For a given frequency $\omega$, we draw thirty realizations of $\mn@boldsymbol{\epsilon}=\cos(\omega\textbf{{t}}+\phi)$ , $\phi$ following a uniform distribution on $[0,2\pi]$. For each of them we compute $R_{\mn@boldsymbol{\epsilon}}$, the least square estimate of the orbital elements $\widehat{\theta}$, the estimate of the noise level $\|\textbf{{y}}-\textbf{{f}}(\widehat{\theta}\|/\sqrt{N}$. The $R_{\mn@boldsymbol{\epsilon}}$ are averaged to have an estimate $\widehat{Q_{\mn@boldsymbol{\epsilon}}}$ of $Q_{\mn@boldsymbol{\epsilon}}$ and the average error on eccentricity $\langle|\widehat{e}-e_{t}|\rangle$. For each type of noise $\mn@boldsymbol{\epsilon}$, we plot $(\widehat{Q_{\mn@boldsymbol{\epsilon}}},\langle|\widehat{e}-e_{t}|\rangle)$, which corresponds to a blue point in Fig. 9. The point obtained with $\mn@boldsymbol{\epsilon}$ being a white Gaussian noise model is represented with a yellow cross. We then bin the values with a constant step in $\log Q_{\mn@boldsymbol{\epsilon}}$ and compute the average error, as well as its standard deviation (purple stair curve). We also compute the mean value of the estimated jitter, and divide it by the value of the jitter estimated for $Q_{\mn@boldsymbol{\epsilon}}=1$. The normalized jitter so obtained is represented in green, with its scale on the right $y$ axis.

Fig. 9 is obtained with the time array of Gl 96 SOPHIE measurements (67 measurements), $e_{t}=0$, an assumed covariance matrix V equal to identity, a period of 40 days and fixed ratio of the norm of the Keplerian signal and the input noise of 10 (which corresponds to a very high SNR = 78), to have as little influence as possible of the noise level. As the noise becomes more correlated with the $\textbf{{M}}_{2}$ space ($Q_{\mn@boldsymbol{\epsilon}}$ increases) it is absorbed in the fit and the RMS of the residual decreases. Indeed in Fig. 9 it is apparent that as $Q_{\mn@boldsymbol{\epsilon}}$ increases the error on eccentricity grows while the estimated level of the noise decreases. The error on $e$ for $e_{t}=0.9$ is consistently higher, which is likely due to local minima at high eccentricity (see section 2.2.3)

One can then wonder how the SNR $S$ affects the bias on eccentricity for correlated noises. Defining the noise “quality factor” as $q:=\sqrt{Q_{\mn@boldsymbol{\epsilon}}(\textbf{{M}}_{2})}$, in the linear approximation used for (11), the uncertainty on $k$ and $h$ becomes $q\sigma_{k}$. With notations of equation (13), the bias at $e_{t}=0$ is therefore approximately

[TABLE]

For a given $q$, as the SNR increases, the bias decreases.

To check if the power at semi period is also correlated with the error on eccentricity if the true eccentricity is high, we perform the same experiment with a value of the eccentricity equal to $0,0.3,0.5,0.7,0.9$. Again, the ratio of the norm of the Keplerian signal and the input noise is fixed to 10. For each eccentricity the experiment is done 20 times with periods randomly drawn with a log-normal law. The results are shown in Fig. 9. We first remark that as $e$ increases, the bias, and therefore the error on eccentricity, decreases. Secondly, from $e=0.5$, $Q_{\mn@boldsymbol{\epsilon}}$ is less relevant to predict the effect of the noise on eccentricity estimates. This result is in accordance with Wittenmyer et al. (2019), who found that planets with e>0.5 are very unlikely to be mistaken for a two circular planet model. As eccentricity increases, the harmonics of order greater than two increase, so that a noise localised at the semi period cannot mimic a higher eccentricity.

The same experiment is performed with a ratio of norm of the signal and the noise of three and an injected circular orbit. We compute $Q_{\mn@boldsymbol{\epsilon}}(\textbf{{M}}_{k})$ with $k=1/2,2/3,1,3/2,2,3,4,5$. The results are plotted on Fig. 10, with a colour code for each $k$. It appears that a strong component of the noise on the harmonic 2, 3, 4 or 3/2 and 5 lead to an increased error on eccentricity (in decreasing order of effect). The noise level is notably underestimated in the $k=2$ case. On the contrary, a strong component of the noise on $k=0$ leads to a smaller error, which is easy to understand. Indeed, since the eccentricity of the injected signal is zero, the noise reinforces the signal.

The results of section 3.4.2 are interpretable with the analysis above. The metric $Q_{\mn@boldsymbol{\epsilon}}$ defined in eq. (27) is computed (here $Q_{\mn@boldsymbol{\epsilon}}(\textbf{{M}}_{2})$) for the noises of section 3.4.2, that are noises generated with nominal uncertainties plus $y_{0},...y_{5}$. For $K=2.5$, 3.5 and 5 m.s-1, we generate 10,000 realisations of these noises, and inject a circular planet at period $P$ with a random phase. We do this simulation for $P=5$ to 100 days per step of 5 days. in each of the $6\times 10,000\times 20$ simulations, we compute the FAP associated with $e=0$ with formula (17), and adopt as a convention that $e=0$ is rejected if the $p$-value is below 0.05. We also compute $R_{\mn@boldsymbol{\epsilon}}(\mathbf{M}_{2})$ as defined in eq. (26). For the 6 noises and 20 periods, we average the values of the $R_{\mn@boldsymbol{\epsilon}}$ to obtain $Q_{\mn@boldsymbol{\epsilon}}$ (27). Fig. 11 shows the proportion of false eccentricity rejected as a function of $Q_{\mn@boldsymbol{\epsilon}}$ for $K=5$ m.s-1. It clearly appears that there is a strong correlation with the power of the noise at the semi period.

Obviously, one can test on a given system if a specific type of noise has a particular impact on a planet with injected parameters, the goal of this section was to identify some generic properties of the noise that lead to spurious inferences.

3.5 Robustness to prior changes

3.5.1 Simulations

In the previous sections we have studied the impact of having a wrong likelihood function. We now turn to the sensitivity of the estimate on the prior probability, here with a numerical experiment. With the formalism of section 3.1, the data are generated with prior and likelihood $p_{t}(\mn@boldsymbol{\theta}),p_{t}(\textbf{{y}}|\mn@boldsymbol{\theta})$, and the analysis is done with $p(\mn@boldsymbol{\theta}),p_{t}(\textbf{{y}}|\mn@boldsymbol{\theta})$. The likelihood is correct, but the prior is not, which corresponds to having an incorrect idea of the population distribution. Note that the prior probability of all parameters have an effect on the eccentricity estimate, but we focus on the prior probability chosen for the eccentricity.

Two distributions are considered. We generate eccentricities according to the distribution Beta(a=0.867, b=3.03) and compute posterior probabilities with a uniform prior. Otherwise the priors are taken as in table 5 and the data are generated as in section 3.4. The results are shown in table 4 as a function of the input amplitude. The errors, measured by MSE and MAE are systematically worse when using the incorrect prior, however with an extra error not exceeding 15%. It seems like an error of $\approx 0.01$ on $e$ could be counted as uncertainty on the prior distribution.

The uncertainty on the prior seems not to be a major concern, at least for the estimation of eccentricity. However, in some cases, one might want to recompute the credible interval with another prior, which is the object of the next section.

3.5.2 Recomputing the posterior without new sampling

The most straightforward way to explore the dependency of the posterior on the prior is to recompute it with another prior distribution. However, this might be lengthy to do it systematically on several systems. We here propose an alternative which consists in multiplying the prior by a constant on a subset of its domain of definition, and to scale it elsewhere. In so doing, the output of the posterior sampler can be used straightforwardly without doing any sampling. In the following we illustrate the process with the prior on eccentricity.

We consider a measurable subset $D$ of $[0,1]$ and its complement $\bar{D}$ in $[0,1]$. For instance an interval $D=[0,e_{0}]$ for some arbitrary $e_{0}\in(0,1]$. Let us denote by $\tilde{\mn@boldsymbol{\theta}}$ the model parameters other than eccentricity, so $\mn@boldsymbol{\theta}=(e,\tilde{\mn@boldsymbol{\theta}})$. We define a new prior $p^{\prime}$ such that for $e\in D$, $p^{\prime}(e,\tilde{\mn@boldsymbol{\theta}})=xp(e,\tilde{\mn@boldsymbol{\theta}})$. To ensure that $\int_{0}^{1}\int p^{\prime}(e,\tilde{\mn@boldsymbol{\theta}}){\rm d}\tilde{\mn@boldsymbol{\theta}}{\rm d}e=1$, we take

[TABLE]

We now want to compute the probability that $e\in C$ for a prior distribution $p^{\prime}(e,\tilde{\mn@boldsymbol{\theta}})$. This one is given by replacing $p$ by $p^{\prime}$ in equation (23). With the notations

[TABLE]

for $E\subset[0,1]$, we compute

[TABLE]

where $p^{\prime}(e|\textbf{{y}})$ is the posterior distribution when the prior is $p^{\prime}(e)$ and $p(\textbf{{y}}|e)$ is the likelihood marginalized on all parameters but eccentricity. When $x=0$, all the prior probability goes to the complementary of $D$ and $\mathrm{Pr}\{e\in D|\textbf{{y}}\}$ goes to 0. If $x=(1-a)/a$, $e$ is certainly in $D$ therefore $\mathrm{Pr}\{e\in D|\textbf{{y}}\}=1$.

The advantage of this calculation is that the integrals can be computed from the posterior samples. Denoting, $|E|$ the number of the MCMC samples that are such that $e\in E\subset[0,1]$ and $N_{\mathrm{tot}}$ the total number of samples, an estimate of $I_{E}$ is $\widehat{I_{E}}=|E|/N_{\mathrm{tot}}$, so that an estimate of (34) is

[TABLE]

The reasoning can be extended straightforwardly to credible regions $D$ and $C$ in the parameter space, and to prior region subdivisions in $D_{1}...D_{q}$ with disjoint $(D_{i})_{i=1..q}$ whose union is the whole parameter space.

Apart from the Markov chain samples, the only quantity needed to use (35) is $a$ as given by (29). This expression might be difficult to compute in general, but in the case where $p(e,\tilde{\mn@boldsymbol{\theta}})=p(e)p(\tilde{\mn@boldsymbol{\theta}})$, $a=\int_{D}p(e){\rm d}e$, which is a one-dimensional integral. Analytic expressions might exist and a Riemann integration is always possible.

Since the integrals to be evaluated from posterior samples are random variables, it must be ensured that they have a controlled uncertainty. When breaking the posterior in many domains $D_{1}...D_{q}$, the procedure outlined may become unreliable if there are not enough independent samples in each $D_{k}$. One can easily compute the effective number of samples in each $D_{k}$, $N_{\mathrm{eff},k}$ by counting how many samples are in that region and dividing by the correlation time-scale. A number of effective samples greater than $n$ gives an accuracy of $\approx 1/\sqrt{n}\times 100\%$ on the probability $\mathrm{Pr}\{e\in D_{k}|\textbf{{y}}\}$. Further investigation is left for future work.

3.6 Model comparison: one eccentric planet or 2:1 mean motion resonance

A system of two planets in 2:1 mean motion resonance can be mistaken for one eccentric planet, and vice versa. We here study the possibility to disentangle the two cases via Bayes factor as a function of the SNR. Two models are considered, $\mathcal{M}_{e}$ and $\mathcal{M}_{1:2}$, respectively an eccentric planet and two circular planets in mean motion resonance:

[TABLE]

Denoting by $g(\tilde{\mn@boldsymbol{\theta}})$ a deterministic model encapsulating other planets, offsets, trends etc.. We also let vary a jitter term $\sigma_{J}$ as in eq. (8). The Bayes factor of the two models is defined as

[TABLE]

where

[TABLE]

with $\mn@boldsymbol{\theta}=(K,k,h,P,M_{0},\tilde{\mn@boldsymbol{\theta}},\sigma_{J})$ or $\mn@boldsymbol{\theta}=(K_{1},P_{1},M_{01},K_{2},P_{2},M_{02},\tilde{\mn@boldsymbol{\theta}},\sigma_{J})$ for $\mathcal{M}=\mathcal{M}_{e}$ and $\mathcal{M}=\mathcal{M}_{2:1}$ respectively. We expect the two models to be distinguishable if the amplitude of the second harmonic of the signal can be resolved (Anglada-Escudé et al., 2010, see eq. 5).

In order to determine at which SNR the Bayes factor allows to disentangle resonant planets and eccentric ones, we perform a numerical experiment. We select a semi-amplitude $K$ and a period $P$, then generate a Keplerian signal with $e=0.25$, random $M_{0}$ and $\omega$. $K$ is chosen on a grid (2,5,8,11,14 m.s-1), to see how the ability to disentangle scenarios evolves with the true SNR.

On the other hand, we generate a two planet system. The outer planet has period $P$ with random phase and semi-amplitude $K_{1}=K$. The inner planet is circular with $K_{2}=K_{1}/4$, so that $K_{2}/K_{1}=e$ of the single planet. The phase of the second planet is chosen uniformly. The simulation is performed in two different settings. In the first one, the period ratio is fixed to $P_{2}=P_{1}/2$. In the second, $1/P_{2}$ is chosen uniformly between $[(1-\alpha)2/P_{1},(1+\alpha)2/P_{1}]$ with $\alpha=0.1$. The rationale behind the choice of $P_{2}$ is that the period ratios of Kepler planets are located within a neighbourhood of 2:1 (Steffen & Hwang, 2015).

The priors chosen to compute the Bayes factor are summarized in table 6. The Bayes factor is computed with the nested sampling algorithm PolyChord (Handley et al., 2015a, b). The performance of the algorithm was checked on the data sets of the Evidence Challenge (Nelson et al., 2018). For each data set, the algorithm is ran at least five times. The $\ln\mathcal{Z}$ estimate is taken as the median of the different runs and the error bars are given by the variance of the empirical median as provided by Kenney & Keeping (1962). The error on the log Bayes factor $\ln$BF = $\ln\mathcal{Z}_{2}-\ln\mathcal{Z}_{1}$ is taken as $(\sigma_{\mathcal{Z}_{1}}^{2}+\sigma_{\mathcal{Z}_{2}}^{2})^{1/2}$.

Fig. 13 and  13 show the results of the simulation, where the measurement time arrays are those of CoRoT-9 (Bonomo et al., 2017b) and Gl 96 (Hobson et al., 2018). The $\ln$ of the Bayes factor is represented, such that the correct model is always at the numerator. The black dashed lines indicate a Bayes factor equal to 150 and 1/150. Bayes factors above 150 and below 1/150 correspond respectively to very strong evidence in favour or against the correct model.

Fig. 13 and  13 show that distinguishing the resonant and eccentric models is possible, especially if the period of the inner planet varies, which is the case in Kepler data. Nonetheless, at low SNR, there are cases of decisive Bayes factor against the correct model. In the low SNR regime, the evidence is dominated by the prior, and the parameter space is hard to explore. Both effects can account for these spurious results.

We point out that the values reported in Fig. 13 and  13 are very sensitive to the prior on $K$. The narrower it is, the “cheaper” it is to add a planet, such that the two planets model is be favoured. The same experiment as above is done with a flat prior on $K$ on $[0,40]$ m/s. In that case, the two planet model is systematically favoured, as shown on Fig. 14.

As a conclusion, it seems good practice to check the influence of the prior on $K$. Secondly, since we expect that as the period ratio of resonant planets not to be exactly one half, if the two planet hypothesis is true then it will be strongly favoured by the Bayes factor. It therefore seems reasonable to consider the eccentric planet as the null hypothesis and not to reject it if there is no strong evidence for the two planet model.

4 Detecting model errors: residual analysis

4.1 Objective

In the previous section, we have seen that adjusting a jitter term is satisfactory in most cases, however we did not envision all possible errors. It is good practice to check if the models considered are plausible at all. One approach to take into account model uncertainty is to consider many models (Jones et al., 2017) and rank them via cross validation, Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC), or even Bayes factor if possible. Alternately, we can test the hypothesis that the residuals are consistent with the model in an absolute sense. This problem is sometimes referred to as the goodness-of-fit problem, and is in general difficult (see Lehmann & Romano, 2005, chap. 14).

We reason as follows: if the set of models is appropriate to describe the data, then the residuals of the best fit must verify certain properties. If they do not, then we reject the hypothesis that there is one acceptable model that explains entirely the data set, among the set of models considered. Furthermore, we would like to obtain hints on the origin of a model misspecification. We expect outliers to change the distribution of the residuals, and astrophysical or instrumental noise to introduce correlations. As a consequence, we consider two types of checks: is the distribution of the residuals approximately Gaussian? and: is there remaining time correlations in the residuals?

4.2 Distribution of the residuals

As in section 2.2.2, we first examine the linear case and show that the results are still helpful in the non linear setting. Let us suppose that we have a linear model $\textbf{{y}}=\textbf{{A}}\mn@boldsymbol{\theta}+\mn@boldsymbol{\epsilon}$ where A is a $N\times p$ matrix and $\mn@boldsymbol{\epsilon}$ is a Gaussian noise of covariance matrix $\textbf{{V}}=:\textbf{{W}}^{-1}$. Let us denote by $\widehat{\textbf{{y}}}$ the least square fit model, and suppose the model ($\textbf{{A}},\textbf{{V}}$) is known. Then the weighted residual

[TABLE]

is a vector of $N$ random variables that are approximately independent, Gaussian of null mean and variance one. To obtain a weighted residual that is a vector of independent Gaussian variables, let us define Q, the matrix such that $\textbf{{J}}=\textbf{{Q}}^{T}(\textbf{{I}}_{N}-\textbf{{W}}^{1/2}\textbf{{A}}^{T}(\textbf{{A}}^{T}\textbf{{W}}\textbf{{A}})^{-1}\textbf{{A}}^{T}\textbf{{W}}^{1/2})\textbf{{Q}}$ is diagonal (it exists). Then the re-weighted residual $\textbf{{r}}_{QW}^{\prime}=\textbf{{Q}}^{T}\textbf{{W}}^{1/2}(\textbf{{y}}-\widehat{\textbf{{y}}})$ has $p$ null components. The $N-p$ others are Gaussian variables of mean 0 and variance 1. In what follows, we denote by $\textbf{{r}}_{QW}$ the vector made of the $N-p$ components of $\textbf{{r}}_{QW}^{\prime}$. These two results are proven in Appendix E.

In practice, A and V are unknown, and we choose models $\textbf{{A}}^{\prime}$ and $\textbf{{V}}^{\prime}$. The two above properties can be used to test if $(\textbf{{A}},\textbf{{V}})=(\textbf{{A}}^{\prime},\textbf{{V}}^{\prime})$ because if so, then the weighted residuals $\textbf{{r}}_{W}$ and $\textbf{{r}}_{QW}$ have a known distribution.

We compute an experimental cumulative distribution function (CDF) of $\textbf{{r}}_{W}$ and $\textbf{{r}}_{QW}$. If our model is correct, then it should be close to the CDF of a Gaussian variable of mean zero and variance one.

4.3 Correlations in the residuals

The test suggested in the previous section is relevant to check the distribution of the residuals without temporal information, and thus is not most adapted to spot correlations. We here adapt the variogram (Matheron, 1963) for unevenly sampled time series, similarly to  Baluev (2013a). The quantity $d(t_{i},t_{j})=\textbf{{r}}_{W}(t_{i})-\textbf{{r}}_{W}(t_{j})$ is plotted as a function of $t_{i}-t_{j}$ for $t_{i}>t_{j}$. If $\textbf{{r}}_{W}$ is indeed independent and Gaussian, $d(t_{i},t_{j})$ should not depend on the time interval.

Secondly, we consider $n$ time bins with constant spacing in $\log t$. For each bin $B$, we compute the sample variance of the $d(t_{i},t_{j})$ such that $t_{i}-t_{j}\in B$. We expect that if there are correlations, these variances should grow as $t_{i}-t_{j}$ increases. To obtain an error bar on the variances, we add an independent Gaussian noise of mean 0 and variance 1 to $r_{W}$ and re-compute the variances for the same time bins. One can alternately add a Gaussian noise of covariance $\mathbf{V}$ to the data, and re-compute the residuals. The error bars are taken as $\pm\sigma$ where $\sigma$ is the standard deviation of the variances estimates per bin.

4.4 Example

Let us now show how it can be used in practice. We take the 214 measurement times of Proxima b Anglada-Escudé et al. (2016). A is made of six columns as defined in Appendix A and fix $\textbf{{x}}_{t}$. We then generate three series of a thousand realisation of $\textbf{{y}}=\textbf{{A}}\textbf{{x}}_{t}+\mn@boldsymbol{\epsilon}$. The covariance matrix of the noise has a kernel $\operatorname{e}^{-|\Delta t|/\tau}$ where $\Delta t$ is the duration between two samples. The three series are generated with a noise time-scale $\tau=0$, 10 and 100 days. For each of the $3\times 1000$ signals generated, we compute the least square fit with the correct matrix A, but with a weight matrix W equal to identity, so our model is entirely correct only in the first case. First, we pick randomly one realization among the 1000 available in each series, and perform the first test whose result is plotted Fig. 16. One clearly sees a pattern: the higher the correlation, the smaller is the difference between residuals. We then consider 6 time bins, and compute the variances of the data and their uncertainties within each bin. The results of these calculations are represented by the stair curves. For correlated noises, the variance increases with the time interval, while it stays compatible with a constant for the white noise. Fig. 16 shows the 1000 empirical CDFs in the three cases.

The plots 16 and 16 are useful indicators of remnant correlations in the residuals and non Gaussianity. However, they do not constitute metrics with known statistical properties. One can potentially test the hypothesis that $\textbf{{r}}_{QW}$ is a realization of such a law with a Kolmogorov-Smirnov test or other metrics such as Anderson & Darling (1954), Shapiro & Wilk (1965), etc. We have tested the Anderson-Darling metric which did not show a high statistical power (see Hara 2017, thesis) and that is not discussed further. The visual inspection, though less quantifiable, seems more accurate.

Study of correlations in RV residuals have already been undertaken for instance by Baluev (2011, 2013b) with a smoothed residual periodogram (Baluev, 2009). The systematic comparison of these statistics is left for future work.

5 Application: 55 Cancri

To illustrate the methods above, the estimation of the eccentricities of the 55 Cnc system is discussed. This system has been extensively studied and many measurements, from several spectrographs, are publicly available: Hamilton (Marcy et al., 2002; Fischer et al., 2008), ELODIE (Naef et al., 2004), HRS (McArthur et al., 2004; Endl et al., 2012), HIRES (Butler et al., 2017). We here focus on the HIRES (Butler et al., 2017) data set, from now on denoted by HIRES17, which illustrates the importance of the noise model choice for reliable eccentricities.

The HIRES17 data set is made of 607 velocity measurements, with on average 4 data points per night. Fig. 17 shows these data with nominal error bars. The inset is a zoom on five points taken within a time interval of 50 minutes at reduced Julian day (RJD) 55351. We first consider the raw data. The posterior distributions of the orbital elements is computed with a model including the five known planets plus a free jitter term as defined in section 3.2, and priors defined as in table 5. We now focus on 55 Cnc f, orbiting at 260.9 days with a minimum mass of $0.1503^{+0.0076}_{-0.0076}$ $M_{J}$ and a period of $259.88^{+0.29}_{-0.29}$ days (Bourrier et al., 2018). In Fig. 18, the posterior distribution of the eccentricity is represented in yellow. The estimate is 0.5 with a 68% credible interval equal to $[0.44,0.58]$. This estimate is very different from the one in Bourrier et al. (2018), $e=0.08^{+0.05}_{-0.0.04}$, where orbital elements estimates are based on all the available radial velocity data.

To check for correlations in the residuals, these are studied with the method of section 4.3. We define $\mathbf{W}=(\mathbf{V}+\widehat{\sigma}_{ML}^{2})^{-1}$, where $\mathbf{V}$ is the nominal covariance matrix, and $\widehat{\sigma}_{ML}^{2}$ is the maximum likelihood estimate of the jitter. The weighted residuals $r_{W}$ as defined in eq. (40) are computed. For all combinations of measurement times $t_{i}>t_{j}$ we represent $d_{ij}:=r_{W}(t_{j})-r_{W}(t_{i})$ as a function of $t_{j}-t_{i}$. We then compute the variance of the $d_{ij}$ such that $t_{j}-t_{i}$ is in a certain time bin. Fifteen such intervals are considered, with a constant length in $\log$ scale. The uncertainty on this variance $\sigma^{2}$ is estimated by bootstrap. The results are shown in Fig. 19, where the blue points represent a couple $(t_{j}-t_{i},d_{ij})$, and the red stair curve represents the variance per bin and its uncertainty. As the time difference grows, the residuals are more and more dispersed, which is indicative of correlated noise .

To account for this correlated noise, we adopt two strategies. First, we include a red noise term in the model of section 3.2. This one is assumed Gaussian with a covariance kernel

[TABLE]

and $\sigma_{R}$ and $\tau$ are included in the posterior sampling with a normal prior on $\ln\tau$ of mean and variance 1. The resulting eccentricity posterior for 55 Cnc f is shown in red on Fig. 18. The posterior median of $\tau$ is 2.3 days, with a 68% confidence interval of $[0.5,9.6]$ days. The second strategy we adopt is to bin the velocity measurements per day. Such a strategy is actually what likely lead to the measurement pattern chosen, as binning data per night averages out stellar oscillations (Dumusque et al., 2011). The posterior so obtained is represented by the white histogram on Fig. 18. In both cases it appears that that the eccentricity is in fact much less constrained. Note that the posterior distribution favours a null eccentricity. We attribute this to the choice of the beta prior, which favours small eccentricities, but also to the fact described in section 2.3: high eccentricity models represent a larger volume of models, and are therefore be penalized by Bayesian estimates.

In section 3.4.3, we argued that the main characteristic of a noise that influences eccentricity estimation is, for a particular time sampling, its power at the semi period of the planet. The results above are indeed interpretable in this term. The quantity $Q_{\mn@boldsymbol{\epsilon}}(\mathbf{M_{2}})$ defined in eq. (27) is computed for noises with a kernel defined in eq. (41). The calculation is made with the nominal time sampling and the one corresponding to a one day binning. Fig. 20 shows the result of these calculations. With the nominal times, the power at 130 days is much more sensitive to correlated noise, even at time-scales as short as 5 minutes, which results in a clear eccentricity overestimation. This example illustrates that due to the sampling, some noises might create spurious eccentricities if not properly accounted for.

As a remark, some measurement dates correspond to the observation of the 55 Cnc e transit. As noted in Bourrier et al. (2018), the effect of the Rossiter-McLaughlin effect has been put in question in López-Morales et al. (2014), and should not exceed 0.5 m.s-1, which is well below HIRES precision and does not impact the analysis above.

6 Conclusion

6.1 Summary and step-by-step method

In section 2, we have seen that at low eccentricity, the bias of the least square estimate can be approximated analytically (eq. (12)). This equation as well as numerical simulations show that the bias is approximately proportional to the uncertainty on eccentricity. We show that an incorrect estimate of the noise level can create spurious global minima at high eccentricity through the Proxima b example, and that at high SNR it is less likely, though not impossible, to miss the global minimum by doing a least square fit initialized at $e=0$. In section 2.3, it appeared that the maximum of the marginal distribution of eccentricity (eq. (20)) is less biased than the maximum likelihood.

In section 3, we explored the reliability of the maximum likelihood and posterior distributions when the noise model is incorrect, along with the sensitivity of the inference to numerical methods. In summary, we recommend the following data analysis method for robust inference of eccentricities.

• •

Computing the posterior distribution of eccentricity leads to correct inference in general (see section 2.3.2), if one includes in the noise model at least a free jitter term (eq. (24)). Without such a term the inference is very likely to be spurious. The mean and median of the marginal posterior distribution of eccentricity, defined in section 2.3, constitute good point estimates (see section 3.4.2). Credible intervals (eq. (23)) allow reliable hypothesis testing, confidence intervals (see section 2.2.4) are an alternative for probing the low eccentricity region especially.

• •

Using a white, Gaussian noise model might however lead to spurious inferences in some cases. Eccentricity estimates are weakly sensitive to non Gaussianity of the noise, but are impacted by its true covariance, especially by the noise power at $P/2$ and to a lesser extent at $P/3$ where $P$ is the period of the planet of interest (see section 3.4.3). Stellar signals might have a modulation at the rotation period, but as shown in section 5, correlated noise combined with the sampling can also create spurious eccentricities. Inappropriate models leave signatures in the residuals, these ones can be studied with techniques presented in section 4.

• •

It is key to check the convergence of the MCMC used for posterior calculation, for instance with the effective number of independent samples (we recommend at least 10,000). An unreliable numerical method increases the bias (see section 3.3).

• •

At low SNR, posteriors are dominated by the prior and the parameter space is harder to explore. This might lead to spurious inferences (e.g. in section 3.6). The influence of the prior can be assessed with the method of section 3.5.2.

• •

The most degenerate case is an inner planet in 2:1 resonance. As discussed in section 3.6, Bayes factors can disentangle those cases. However, the value of the Bayes factor strongly depends on the prior used for the semi amplitude. We suggest to select a $\log$ prior on semi amplitude and to consider the eccentric scenario as the null hypothesis.

• •

In multiplanetary systems, checking the system stability might also help ruling out some values of the eccentricity (e.g. Hébrard et al., 2016; Delisle et al., 2018).

• •

As more measurements are obtained, least square and Bayesian estimates get closer, so that eq. (13) can also be used to approximate the number of measurement $N(\epsilon_{e})$ needed to have an average bias on eccentricity $\epsilon_{e}$ for a planet with semi-amplitude $K$, $p$ parameters fitted in total and a measurement root mean square RMS,

[TABLE]

This formula is for $\epsilon_{e}\leqslant 0.05$ and subtends that there is a good phase coverage.

This procedure should be familiar to observers, since it is mainly a formalized and tested version of common practices. It is applicable to other purposes, especially the methods of section 4, since checking the model validity improves the inference robustness in general.

6.2 Perspectives

The main point of the present work is that modelling errors might have systematic impacts on the estimates of orbital elements, and thus on our understanding of planetary systems. It is yet to be determined in which extent this has been the case in past studies.

In section 4, we presented tools to measure the absolute adequacy between a model and the data, which is complementary to comparing models to one another. There remain many such adequacy metrics to explore. These ones could prove useful in the context of exoplanets characterization but also exoplanets detection.

Acknowledgements

N. Hara and J.-B. Delisle acknowledge the financial support of the National Centre for Competence in Research PlanetS of the Swiss National Science Foundation (SNSF). N. Hara thanks J. J. Zanazzi for his interesting inputs. We thank the anonymous referee for his/her insightful suggestions.

Appendix A First order approximation

In this section, we developp the Keplerian model to first order in eccentricity to obtain an analytical expression of the bias. Within this approximation, the distribution of the least square fit knowing that $e=0$ is given in Lucy & Sweeney (1971). This section extends their formula to small $e$, and takes into account the number of fitted parameters. First, we develop (4) to order one in $e$, obtaining

[TABLE]

where $\lambda=nt+\omega+M_{0}=\lambda_{0}+2\pi t/P$ is the mean longitude, $\lambda_{0}$ being its value at $t=0$. Denoting by $n=2\pi/P$ the mean motion, the above expression can be re-written

[TABLE]

where $A=K\cos\lambda_{0}$, $B=-K\sin\lambda_{0}$, $C=Ke\cos(2\lambda_{0}-\omega)$, $D=-Ke\sin(2\lambda_{0}-\omega)$. When other parameters are fitted, the uncertainties on $A,B,C,D$ increases as well. To quantify this effect, we consider the problem of fitting the period and a constant.

[TABLE]

which in a matrix form gives

[TABLE]

Let us assume that the observations are $\textbf{{y}}(\mathit{t})=\textbf{{M}}(P)\mathit{x}_{t}+\mathit{\epsilon}$, where $\mathit{\epsilon}$ is a Gaussian noise, independent and identically distributed with variance $\sigma^{2}$. The least square estimate of x is $\widehat{\textbf{{x}}}=(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}\textbf{{y}}$, and the estimate of eccentricity is

[TABLE]

where $K_{t}$ is the true semi-amplitude. By change of random variable we can obtain the law followed by $\widehat{e}$. If we assume that $N$ is large enough then the columns of $\textbf{{M}}(P)$ are approximately orthogonal, the components of $\widehat{x}$ are independent Gaussian variables. Since the modulus of a sum of independent Gaussian variables follows a Rice distribution,

[TABLE]

where $I_{0}$ is a modified Bessel function of first kind, $S=K_{t}/\sigma$ is the SNR, where $\sigma$ is the standard deviation of $\widehat{A},\widehat{B},\widehat{C}$ and $\widehat{D}$. If $K$ is sufficiently large, $W$ is close to 1 and $g(u)$ gives a good approximation of the law followed by the eccentricity fitted. Within this approximation, one can obtain analytical formula for the bias $b$ of the eccentricity that only depends on the true eccentricity and the SNR,

[TABLE]

where $L_{1/2}$ is the Laguerre polynomial of order $1/2$. In case $K_{t}$ is small, one must use the formula for the law followed by the quotient of two random variables:

[TABLE]

but no simple analytical expression was found.

When fitting model (46) to $\textbf{{y}}(\mathit{t})$, the estimate $\widehat{\mathit{\theta}}$ have a covariance matrix $\bm{\Sigma}^{-1}$ where $\bm{\Sigma}=\sigma^{2}(\textbf{{M}}(P)^{T}\textbf{{M}}(P))$ (this is a classical statistical result, see for example Pelat (2013)). The variances of the components of $\widehat{\mathit{x}}$ are given by the diagonal elements of $\bm{\Sigma}^{-1}$. Their approximate calculation is the object of the next section.

A.1 Average error

First we consider the estimation of the error on $A,B,C,D$ when averaging over the mean motion $n$. At little cost, we can generalize our claim to the fitting of model (45) plus fitting other linearised Keplerian model. This approximately corresponds to fitting a multi-planetary system starting closely from the correct local minimum of $\chi^{2}$. Again, the model can be written as a linear one, $\textbf{{y}}=\textbf{{M}}\textbf{{x}}$ but where M has $p=6+5k$ columns, $k$ being the number of additional planets. To facilitate the discussion, we normalize the columns of M. To have an expression of the model of the form (46), we have multiplied the $k$th component of $\theta$ by the norm of the $k$th column of $M$. The variances of these new model parameters are still given by the diagonal elements of $\sigma^{2}\bm{\Sigma}^{-1}$ where $\bm{\Sigma}=(\textbf{{M}}^{T}\textbf{{M}})$, but now $\bm{\Sigma}$ has only ones on its diagonal.

Calculating precisely the uncertainty on $A,B,C,D$ averaged over $n$ and the phase of the signal as a function of the instant of observations $t$ is complex since it requires the inversion of $\sigma$ which is a $6+5k\times 6+5k$ matrix. Instead, we use an approximation that grasps the effect we want to estimate: how the uncertainty worsens as more parameters are added to the model. We consider that the elements of M are drawn from independent Gaussian laws that have a variance $1/N$. To avoid confusion with the true model, the so defined random matrix is denoted by $\tilde{\textbf{{M}}}$ and its covariance matrix by $\tilde{\bm{\Sigma}}$.

This approximation seems to be rough at first but turns out to be surprisingly accurate as a lower bound in practice. A few arguments to justify that it is a reasonable guess are listed below.

• •

The variances of the entries were chosen such that the expectancy of a squared norm of a column is one, which is the value of $\bm{\Sigma}$ diagonal elements.

• •

The columns are cosines and sines, which are approximately orthogonal, and in the Gaussian case decorrelation implies independence. Furthermore, the average of the spectral window is equal to the expected value of a correlation between two random Gaussian variables.

• •

The normed vectors $\cos\nu\mathit{t}$ and $\sin\nu\mathit{t}$ are approximately distributed uniformly on the sphere of $\mathbb{R}^{N}$ when $\nu$ is distributed uniformly between 0 and $2\pi/T_{\mathrm{obs}}$.

The expected value of the variance of any parameter is the expected value of any diagonal element of $\bm{\Sigma}^{-1}$, since all the columns of $\tilde{\textbf{{M}}}$ follow the same law. To tackle that problem, we rewrite $\tilde{\bm{\Sigma}}$ as

[TABLE]

Where $\bm{\Sigma}_{11}$ is $\tilde{\bm{\Sigma}}$ element at first row and first column and $\bm{\Sigma}_{1}$ is a column vector with $N-1$ entries. We now have

[TABLE]

By Jensen inequality (Jensen, 1906), since $x\rightarrow 1/x$ is convex,

[TABLE]

Now since for two independent variables $X$ and $Y$, $\mathbb{E}\{XY\}=\mathbb{E}\{X\}\mathbb{E}\{Y\}$,

[TABLE]

As by construction $\mathbb{E}\left\{\bm{\Sigma}_{11}\right\}=1$, we finally obtain

[TABLE]

where the inequality follows again from Jensen’s inequality applied to matrix inversion. Finally, the standard deviation on $I=A,B,C,D$ is

[TABLE]

With the approximation $\|\cos\nu\textbf{{t}}\|\approx\|\sin\nu\textbf{{t}}\|\approx\sqrt{N/2}$, the errors on $k=C/\sqrt{A^{2}+B^{2}}$ and $k=D/\sqrt{A^{2}+B^{2}}$ then verify

[TABLE]

As $k$ and $h$ approximately follow a Gaussian law, $e=\sqrt{k^{2}+h^{2}}$ follows a Rice distribution, whose mean is given by

[TABLE]

$L_{1/2}$ being the Laguerre polynomial of degree 1/2. The relevance of formula (54) is checked on numerical examples next section. As we shall see, the lower bound is tight when $p$ does not exceeds $\approx N/2$.

Let us finally stress that formula (54) approximates the bias averaged on the mean motion, that is the frequency of the orbit, not the period. Averaging on the period would give more weight to the bias at low frequencies, which is high, and would therefore lead to a greater average value of the bias.

A.2 Precision and accuracy

Eq. (3) expresses the mean squared error (MSE) as a function of the bias and the standard deviation. Assuming the model is correct, the MSE is an accuracy metric (dispersion about the true value), while the variance captures the precision of the estimate (dispersion of the estimate about its mean value). Shen & Turner (2008); Zakamska et al. (2011) have shown that the estimates accuracy degrades – all other parameters being fixed – as the SNR decreases, as the period of the planet is longer, and as phase coverage degrades.

In this appendix, we argue that theses effects can be seen as degrading the precision of the eccentricity estimates. More precisely, we show that the bias is proportional to the standard deviation of the estimate, so that the MSE is also proportional to the standard deviation.

We proceed with a numerical simulation. We consider the 74 measurement times of HD 69830 (Lovis et al., 2006), spanning on 800 days, as they are spaced in a typical manner. We then inject a simulated planet in circular orbit and a white, Gaussian noise of standard deviation 1 m.s-1. In all the following simulations, the phase is uniformly random. By default, the semi-amplitude of the the planet is $K$ = 3 m.s-1 and the period is 31.56 days period (like HD 69830 c).

Data points are taken off two by two, from 74 to 14. 2. 2.

The semi amplitude of the planet varies from 0.5 to 6.3 m.s-1 by a step of 0.2 m.s-1. 3. 3.

The period is drawn from a log-normal law, where $\log_{10}P\sim G(1,1)$. Thirty different periods are drawn. 4. 4.

The phase coverage is degraded. We consider the 25 times $t_{k}=k\times P$ where $P=31.56$ is the planet period. For each $t_{k}$, three epochs of measurements are drawn uniformly between $t_{k}$ and $t_{k}+\Delta t$. We choose thirty different lengths for $\Delta t$, equispaced from $P/4$ to $P$. The rationale is to generate observations more or less localised around the same time when folded in phase at $P$.

Each of the four simulations is made with thirty different parametrizations. For each of the $30\times 4$ of these, we generate 500 realisations of white noise and report the average value of eccentricity (the bias) and the standard deviation of the estimate. The results are reported in Fig. 21, where it appears that the bias is proportional to the standard deviation. The points obtained with simulations 1 to 4 described above correspond to the blue, red, yellow and purple points on the graph. The analytical approximation (13) suggests that the bias should be proportional to the standard deviation of the estimate with a factor $\sqrt{\pi/(4-\pi)}$. The black line, which represents $y=\sqrt{\pi/(4-\pi)}x$ is in close agreement with the scatter observed.

The $\sigma_{e}$ reported in Fig. 21 is computed as the standard deviation of the estimate. Note that the bias is also proportional to the uncertainty on eccentricity computed from the correlation matrix, obtained from the least square fit.

Appendix B Local $\chi^{2}$ minima at high eccentricities

As shown in Baluev (2015), the number of local minima increases significantly in the high eccentricities region. These minima might lead a local minimisation algorithm or a Monte Carlo Markov Chain (MCMC) to be stuck in the wrong region of the parameter space. We here aim at quantifying and understanding this feature.

In section 2.2.2, we defined a notion of SNR (eq. (11)). Interestingly enough, one can adapt this notion to predict the number of local minima at high eccentricity, provided the value of the SNR is checked to be reliable.

B.1 Number of local minima per SNR

We explore through simulations how many local minima should be expected, and how frequently the local minimum given by a non-linear fit starting at a circular orbit is not the global minimum. Our simulation is structured as follows. We consider the measurement times of CoRoT-9 (Bonomo et al., 2017b), Gl96 (Hobson et al., 2018), of HD69830 (Lovis et al., 2006), HD40307 (Mayor et al., 2009), which have respectively 28, 67, 74 and 129 observations and generate a planet with uniform distribution in $e$, $\omega$, $M_{0}$ and $K$ (the latter on [0,3]). The period is drawn from a log normal distribution ($\log_{10}P\sim G(1.5,1)$). The uncertainties on the measurements are taken as the nominal ones that are normalised to obtain a mean variance of 1. We then obtain the fitted uncertainty $\sigma_{\mathrm{fit}}^{2}$.

For each realisation of the true planet and the noise, we perform first a local minimisation, initialized at a circular orbit. We then compute three quantities at this local minimum: the condition number of the Fisher matrix F, as well the SNR, similarly to equation (11),

[TABLE]

where $K_{t}$ and $K_{\mathrm{fit}}$ are respectively the true and fitted values of the semi-amplitude, $N$ is the number of measurements, $p$ the number of fitted parameters. $\sigma_{\mathrm{fit}}$ is defined as follows. We adjust a term $\sigma_{J}^{2}$ so that the reduced $\chi^{2}$ is equal to one and take $\sigma_{\mathrm{fit}}$ as the mean of $(\sigma_{J}^{2}+\sigma_{k}^{2})_{k=1..N}$. The rationale of taking this definition instead of the true SNR is to have a quantity that does not require to know the true orbital elements and can be computed on a real data set.

Secondly, we perform a Keplerian periodogram on a closely spaced grid of $e\in[0,1]$, $\omega\in[0,2\pi)$ and $n\in[n_{t}-1/T_{\mathrm{obs}},n_{t}+1/T_{\mathrm{obs}}]$, where $n_{t}$ is the true mean motion. We compute the value of the eccentricity where the global minimum of $\chi^{2}$ is attained, as well as the number of local minima. For instance, in Fig. 3, the red curve displays four local minima at $e=0.17,0.75,0.92$ and $1$, and the global minimum is attained at $e=0.92$. In that case, a local minimisation would not give the global minimum.

For both definitions of the SNR, we compute the number of system that have a SNR between $5k$ and $5(k+1)$, $k=0..7$ and greater than 35. In each bin, we compute the number of systems with 1,2,3,4,5,6 or 7 local minima. The results are represented in Fig. 23 and 23 for $S_{t}$ (eq. (57)) resp. $S_{\mathrm{fit}}$ (eq. (58)) by the blue histograms. We then compute in each bin the fraction of system where the global minimum is not the local minimum closest to 0, which means that most likely, a local minimisation does not yield the global minimum (red stair curve, see the right $y$ axis scale). For instance, in Fig. 23, we see that out of the 2500 systems simulated, 472 had a St between 0 and 5, 176 of which had only one local minimum, 203 had two, 84 had 3, 9 had 4 and none had more. Among those 472 systems, 34% of them had a local minimum that is not the global minimum.

Note that in Fig. 23, the systems with $S_{\mathrm{fit}}>35$ present the highest proportion of missed global minima. This is due to the fact that no systems with very high SNR where generated. As a consequence, all the very high values of $S_{\mathrm{fit}}$ result from datasets with very low $S_{t}$ where $K$ was very overestimated. To obtain a more reliable diagnostic, we need to determine if a fitted SNR can be trusted. The criterion we used is to select only the data sets where $c_{F}<10^{7}$. The figure 24 shows the number of local minima per bin of $S_{\mathrm{fit}}$ so obtained. The predominance of systems with two local minima is due to the fact that at very high eccentricity, there is in general a decrease of $\chi^{2}$. As a consequence, the second local minima is attained at the maximum eccentricity used for the calculation of the global periodogram. For instance in the case of Proxima b (Fig. 3), there is such a minimum at $e=0.999$. Note if we do the same analysis per data set (CoRoT-9, Gl96, HD69830, HD40307), we obtain very similar figures, which shows that indeed the SNR is a reliable metric for the number of local minima. As a conclusion, we expect that the exploration of the parameter space to be more difficult at low SNR, since there are more local minima to expect.

B.2 Interpretation

We now give a geometrical interpretation of why the high eccentricity region is prone to having local minima.

Finding the best fitting model amounts to finding the model closest to the observation in a geometrical sense. We consider the figure drawn in $\mathbb{R}^{N}$ by all the models that have an eccentricity $e$ and a period $P$, $\mathcal{M}_{e,P}$. This figure might explore more or less dimensions. For instance, if it is close to a plane, it is nearly confined to a two-dimensional space. Otherwise, exploring many dimensions traduces a “rough” surface, which increases the chances of finding a local minimum of distance to the data.

By a procedure based on singular value decomposition, detailed in Hara (2017) section 4.3, we obtain an approximate number of dimensions explored by $\mathcal{M}_{e,P}$ as a function of $e$. We here provide a brief description of our methodology. The measured velocity can be expressed as a linear combination of the velocity components in the orbital frame $\dot{X},\dot{Y}$. With these variables, equation (4) becomes

[TABLE]

The components $\dot{\textbf{{X}}}(P,e,\omega)=(\dot{X}(t_{k},P,e,\omega))_{k=1..N}$, $\dot{\textbf{{Y}}}(P,e,\omega)=(\dot{Y}(t_{k},P,e,\omega))_{k=1..N}$ are computed for a grid of $\omega$, $(\omega_{k})_{k=1..n}$. Those vectors are concatenated to form a matrix $M(e,P)=[\dot{\textbf{{X}}}(P,e,\omega_{1})..\dot{\textbf{{X}}}(P,e,\omega_{n})\dot{\textbf{{Y}}}(P,e,\omega_{1})..\dot{\textbf{{Y}}}(P,e,\omega_{n})]$, whose columns are normalized to obtain $\tilde{M}(e,P)$. We compute the number of singular values of $\tilde{M}(e,P)$ that are greater than a tenth of the maximum singular value. This constitutes a proxy for the number of dimensions explored by models with $e,P$ fixed.

As an example, we perform this calculation on the 214 measurement times of GJ 876 (Correia et al., 2010). The number of explored dimensions are shown in table 2. As eccentricity increases, the models explore more dimensions. Since the models at high eccentricity occupy a very large volume in many dimensions, there will often be at least one high eccentricity model closely fitting the data.

Appendix C Frequentist methodologies

C.1 Presentation

Testing possible eccentricities can also be done in a frequentist framework. This one offers confidence intervals, which are not as easy to interpret as Bayesian credible intervals but have the advantage of being quicker to compute. Furthermore, the associated algorithms have clearer convergence tests.

So far frequentist inferences for eccentricities have been done in several ways. (Lucy & Sweeney, 1971) and Husnoo et al. (2012) respectively used $p$-values and Bayesian Information Criterion (BIC) to test the hypothesis that eccentricity is non zero. More precisely, Lucy & Sweeney (1971) compute the probability distribution of the eccentricity estimate under the hypothesis that the eccentricity is null and find a Rayleigh distribution whose variance depends on the SNR (which we also obtain as a special case of our analysis section 2.2.2 with $p=0$). For a given measured eccentricity $\widehat{e}$, they measure the probability that the Rayleigh distribution is higher than $\widehat{e}$ and report an eccentric orbit if this probability is lower than a certain threshold.

Husnoo et al. (2012) computes

[TABLE]

where $\chi^{2}_{\mathrm{min}}$ is the minimum $\chi^{2}$ obtained when minimising the distance between the data and model $\mathcal{M}$, $p$ is the number of degrees of freedom of $\mathcal{M}$ (three for a sine model and five for a Keplerian one) and $|\textbf{{V}}|$ is the determinant of the correlation matrix. The orbit is said to be eccentric if $\mathrm{BIC}(\mathcal{M}_{\mathrm{ecc}})\geqslant\mathrm{BIC}(\mathcal{M}_{\mathrm{circ}})$ where $\mathcal{M}_{\mathrm{ecc}}$ and $\mathcal{M}_{\mathrm{circ}}$ are respectively eccentric and circular models.

Though reasonable, these techniques can be improved. First, they both consider the alternative $e$ is zero or non-zero, and do not allow to test if a given value of eccentricity is compatible with the data or not. Secondly, the analytical approximation of the eccentricity distribution is not always accurate. Finally, the Bayesian information criterion (60) gives equal weight to all parameters, only their number $k$ appears. This approximation is valid in the limit of a large number of observations.

C.2 New methods

Our aim is to overcome as much as possible these limitations. It turns out that the procedure to construct confidence intervals outlined in Casella & Berger (2001), chapter 9, allows us to test the hypothesis that the true eccentricity is equal to a certain value $e$ for all $e$. The idea is to reject the hypothesis that eccentricity is equal to $e$ if all models with eccentricity $e$ have a likelihood lower than a fraction of the maximum likelihood. The following criterion is computed in Appendix C.4. We reject the hypothesis that the eccentricity has a certain value $e$ with a confidence level $\alpha$ if

[TABLE]

where $\Theta_{e}$ is the set of parameters that have all eccentricity $e$, $f(\textbf{{y}}|\mn@boldsymbol{\theta})$ is the likelihood, $F_{\chi^{2}_{\rho}}^{-1}$ is the inverse cumulative distribution function of a $\chi^{2}$ law with $\rho$ degrees of freedom, $S^{\prime}=(\sigma/K_{t})\sqrt{2/N}$ and $L_{\frac{1}{2}}$ is the Laguerre polynomial of order $1/2$. The quantity (61) is simply the ratio of the maximum likelihood obtained by restriction to the models with fixed eccentricity divided by the maximum likelihood on all models. The condition states that if all models that have eccentricity $e$ have too low a likelihood, then $e$ is rejected. The following equations give the value of that threshold, which is obtained by calculating the law followed by the random variable $LR$ under the hypothesis that the true eccentricity is $e$ ($LR|(e_{t}=e)$). It is in fact easier to compute the law followed by the logarithm of $LR$, to obtain a $\chi^{2}$ law whose degree depends on a definition of the SNR $S$ and on the eccentricity under study, but is always smaller than 2. Our computations, detailed in Appendix C.4, also make use of simplifying assumptions, but these are checked to give satisfactory results on simulated signals.

One of the problems of that expression is that it depends on the true value of the semi-amplitude, $K_{t}$, which is unknown. There are two ways to circumvent this issue: either by assuming that $\rho=2$ for all $e$, which is the maximum value $\rho$ can take, of $K_{t}$, or by approximating $K_{t}$ by the semi-amplitude of a circular orbit fitted at the period of the signal. The first option can be used to obtain conservative intervals to ensure that $e$ is non zero. The second one gives a more realistic criterion to reject an eccentricity if no extra care is needed. Let us note that $\rho=2$ is obtained for $e=0$. This has a simple interpretation: the model can be approximated by a linear one in $k=e\cos\omega$ and $h=e\sin\omega$. When $e=0$, both $k$ and $h$ are set to zero, which blocks two degrees of freedom. Denoting by $\textbf{{y}}_{e}$ the model with fixed eccentricity $e$ that has maximum likelihood and $\textbf{{y}}^{\star}$ the model with maximum likelihood, all parameters free,

[TABLE]

behaves then as a $\chi^{2}$ law with two degrees of freedom.

To offer a point of comparison of the confidence interval computed equation (63), two other metrics are tested. First, we simply consider

[TABLE]

which is basically equation (64), normalized by $\|\textbf{{W}}(\textbf{{y}}-\textbf{{y}}^{\star})\|^{2}$ so that it depends less on the noise level assumption. The quantity (65) is assumed to follow a $F$ distribution with $\rho$ and $N-p$ degrees of freedom.

Secondly, we generalize the test suggested by Lucy & Sweeney. Let us denote by $e^{\star}$ the estimate of eccentricity obtained by maximum likelihood when all parameters are free. For eccentricity $e$, we fit a Keplerian model that has an eccentricity fixed at $e$. We then compute the probability

[TABLE]

that is the probability that an eccentricity estimated by maximum likelihood $\widehat{e}$ is at least as far from its assumed value $e$ than the distance between $e$ and the best fit actually observed, assuming the noise model is Gaussian with known covariance matrix V. We also assume that $e$ follows a Rice distribution as in appendix A. A Rice distribution can be seen as the modulus of a vector with two independent Gaussian variables that have the same variance, $X\sim G(a,\eta^{2})$ and $Y\sim G(b,\eta^{2})$ where $k$ and $h$ are the means of these variables. To specify the distribution, we need therefore two scalars: the variance of both random variables $\eta^{2})$ and the modulus of the mean of these two variables, $r=\sqrt{a^{2}+b^{2}}$. Here $a=k$ and $b=h$, so $r=e$. Then $\eta^{2}$ is the variance of the estimates of $k$ or $h$, which under the hypotheses of section 2.2.2 have the same variance $\eta^{2}=(\sigma_{\mathrm{RV}}^{2}/K_{e}^{2})(\pi/(N-p))$. Then the quantity (66) can easily be evaluated by the cumulative distribution function of the Rice distribution, which is a Marcum $Q$-function.

Computing (64), (65) or (66) necessitates to compute the minimum distance between the observations and a model with fixed eccentricity. To do so, we exploit the fact that Keplerian models are partly linear, $\textbf{{y}}(\textbf{{t}},\mn@boldsymbol{\theta})=A\dot{\textbf{{X}}}(P,e,\omega)+B\dot{\textbf{{Y}}}(P,e,\omega)+C$ where $\dot{\textbf{{X}}}$ and $\dot{\textbf{{Y}}}$ are the components of the velocity on the orbital plane. For each couple $e,\omega$, we can minimize $\|\textbf{{y}}-\textbf{{y}}(\textbf{{t}},\mn@boldsymbol{\theta})\|$ on $A,B,C$ and $P$, which are respectively three linear parameters and one non-linear parameter. Such a problem is fast to solve with, for instance, a Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963). If the period is already known (which is supposed here), obtaining an array of $\chi^{2}$ on a fine grid of $e$ and $\omega$ (60 values each) takes only up to one minute. Let us finally note that the idea of restricting the global $\chi^{2}$ minimization to a grid of non linear parameters is not new (Hartkopf et al., 1989; Lucy, 2014). There are even further resemblances of our interval calculation with Lucy (2014), where confidence intervals on orbital parameters are computed in a similar way. However,  Lucy (2014) uses a degree of freedom $\rho=1$ for all parameters. This is correct only if the model is linear in all the parameters or approximately linear in the vicinity of the best fit and unimodal.

C.3 Tests

The formula (62) and (63) have been derived with simplifying assumptions. To test and compare them to other options, we proceed as follows. We define the acceptable interval as the set of $e$ where $LR_{e}\leqslant\exp(-0.5F_{\chi_{\rho}^{2}}(1-\alpha))$.

We generate a population of exoplanets according to a certain prior density on the orbital elements $p(K,e,P,M_{0},\omega)$. The measurement times are taken from existing data sets. The noise generated according to a Gaussian density of covariance matrix $V$. 2. 2.

For each system, we compute the set of eccentricity that are not rejected, we check that the true eccentricity belongs to this set and compute the measure of its complement in $[0,1]$, that is the measure of the set of rejected eccentricities. 3. 3.

The results are summarized in two plots. First, the fraction of cases where the true eccentricity is not in the acceptable interval as a function of $\alpha$. Second, the curve drawn when $\alpha$ goes from 0 to 1 by a point whose ordinate is measure of the complement of the set of acceptable eccentricity and whose abscissa is the fraction of cases where the true eccentricity is not in the acceptable interval.

Such tests were carried out with the following inputs: the measurement times are those of CoRoT-9 (Bonomo et al., 2017b). The angles $\omega$ and $M_{0}$ are chosen uniformly, $e$ follows an uniform distribution. In Fig. 25 a) and b)we plot the result of the experiment for a period fixed at 95 days the semi-amplitude is fixed to $K=3.5\sigma$ where $\sigma$ is the RMS of the errors. These are the parameters of the detected Jupiter in the system. In Fig. 25 c) and d), we let the period vary uniformly in $\log P$ and compute the same quantities.

Plots a) and c) of Fig. 25 are labelled “ROC - like” curve as a reference to Receiver-Operator Characteristic. These ones are defined when the data are used to decide between two hypothesis. The ROC curve represents the fraction of false positives as a function of the fraction of false negatives for a given decision rule. We adapt this notion to our case, where there is an infinity of hypotheses (each $e$ in $[0,1]$ is a hypothesis). For a given rate of true eccentricity that is not in the acceptable interval (false negatives), the $y$ axis gives the precision on the estimate. The more eccentricities are rejected, the more precise the estimate. The closer such a curve is to the upper left corner the better: regardless of the value of $\alpha$, the fraction of true $e$ rejected is zero (no false negatives) and almost all other eccentricities are rejected: the estimation is very precise.

Interestingly enough, the ROC curve (left) is very similar for all the metric considered with a slight advantage for the $F$-ratio and the likelihood ratio tests (formula (64) and (65)), which have a better precision (more eccentricities rejected) when the fraction of true $e$ rejected is low. We now need to set the level of true $e$ rejected. As expected, the curve obtained for $\rho=2$ gives an overestimated error rate for a given $\alpha$. For the three other tests, the correspondance seems appropriate. Overall, the $F$-ratio and the likelihood ratio tests seem to perform best.

One advantage of the frequentist method is that it relies only on local minimization algorithms, therefore it is fast and convergence is ensured. We have also shown that the parameter $\alpha$ (see eq. (62)) allows us to directly control the confidence intervals meaning. Let us mention that we observed some peculiar behaviour of the estimates for some periods where the matrix of the linearized model is ill-conditioned, that we wish to investigate into more depth in future studies. In those cases, the hypotheses allowing to compute formula (63) are not verified and Bayesian analysis or more sophisticated formula would be required. On the other hand, it seems unlikely that someone would want to prove a non zero eccentricity of a planet particularly poorly sampled.

C.4 Confidence interval calculation

In this section, we outline the calculation of the confidence intervals for eccentricity. Such an interval is constructed as a set of eccentricities that are not rejected by a hypothesis test. We choose the likelihood ratio test:

[TABLE]

where y denotes the actual observations, $f(\textbf{{y}}|\mn@boldsymbol{\theta})$ denotes the likelihood, $\Theta_{e}$ is the set of parameters that have eccentricity $e$, and $\beta$ is a constant which will be made explicit later. Our aim is to compute the distribution of $R$ under the assumption that the random variable giving the observations is $\mathit{Y}=\textbf{{y}}_{t}+\mathit{\epsilon}$, $\mathit{\epsilon}$ being a Gaussian noise. We further assume the noise is independent and identically distributed, the condition translates to

[TABLE]

where $\mn@boldsymbol{\theta}_{e}=\arg\min\limits_{\mn@boldsymbol{\theta}\in\Theta_{e}}\|\textbf{{y}}-\textbf{{y}}(\mn@boldsymbol{\theta})\|^{2}$, $\sigma^{2}$ is the variance of the observations and $\textbf{{y}}^{\star}$ is the global minimum. We now compute the law followed by $D$, so that we can select a $\beta$ that corresponds to a false alarm probability. Since $D$ is defined implicitly, the calculation of its distribution is difficult. We make two simplifying assumptions that allow us to obtain an analytical expression. The expression will then be tested on real cases through numerical simulations.

Let us first consider the linear approximation $\textbf{{y}}=\textbf{{M}}\textbf{{x}}_{t}+\mn@boldsymbol{\epsilon}$ where M is defined as in (45) and (46). We further suppose that the columns of M are orthonormal. Since the columns are originally of the form $\cos n\textbf{{t}},\sin n\textbf{{t}},\cos 2n\textbf{{t}},\sin 2n\textbf{{t}}$, they must be multiplied by $\sqrt{2/N}$ and the amplitude of the signal is no $K_{t}$ but $K_{t}\sqrt{N/2}$. We look for the solution $\hat{\mathit{\theta}}_{e}$ defined as

[TABLE]

Thanks to the Lagrange multipliers theorem, $\hat{\mathit{\theta}}_{e}$ satisfies the conditions

[TABLE]

[TABLE]

with

[TABLE]

The condition $\partial L/\partial\mathit{x}=\mathit{0}$ leads to

[TABLE]

Since the columns of M are orthonormal, ${\textbf{{M}}}^{T}\textbf{{M}}$ is the identity, thus

[TABLE]

and $x_{j}=u_{j},\ \forall j\geq 5$, where we have defined $u_{i}={\mathit{M}_{i}}^{T}\textbf{{y}}$, $\mathit{M}_{i}$ being the $i$-th column of M. The first four components of $\mathit{x}$ are also constrained by $\partial L/\partial\lambda=0$. Let $U=u_{1}^{2}+u_{2}^{2}$ and $V=u_{3}^{2}+u_{4}^{2}$. We get

[TABLE]

or, equivalently,

[TABLE]

whose solutions are

[TABLE]

For the solution $\hat{\mathit{\theta}}_{e}$ to actually be a minimum of $L$, all its eigenvalues must be positive, i.e., $\lambda$ must verify $-1<\lambda<1/e^{2}$. Only $\lambda_{-}$ fulfils this criterion, thus

[TABLE]

and

[TABLE]

with $e_{0}^{4}=e^{2}V/U$. After a few calculation, we show that

[TABLE]

Let $x=e\sqrt{u_{1}^{2}+u_{2}^{2}}/K_{t}$ and $y=\sqrt{u_{3}^{2}+u_{4}^{2}}/K_{t}$. These two random variables follow Rice distributions with parameters

[TABLE]

An expansion of the product term shows that $D$ behaves approximately as a weighted sum of variables following a $\chi^{2}$ distribution. We can then use the Welch-Satterthwaite approximation (Satterthwaite, 1946; Welch, 1947): $D$ approximately follows a $\chi^{2}$ distribution whose number of degrees of freedom $\rho$ is given by $\mathbb{E}\{D\}$. In the following, we denote by $S^{\prime}=\frac{K_{t}}{\sigma}\sqrt{\frac{N}{2}}$ the SNR. The expected value of $D$ is

[TABLE]

With $\rho=\mathbb{E}\{D\}/\sigma^{2}$, we get

[TABLE]

To obtain a confidence level $\alpha$, then we need to take $-2\ln\beta=F_{\chi^{2}_{\rho}}^{-1}(1-\alpha)$ where $F_{\chi^{2}_{\rho}}^{-1}$ is the inverse cumulative distribution function of a $\chi^{2}$ distribution with $\rho$ degrees of freedom. Conversely, it is possible to convert a measured $D$ to a probability simply by computing $\alpha_{e}=1-F_{\chi^{2}_{\rho}}(D)$. The hypothesis $e_{t}=e$ is rejected if $\alpha_{e}$ is below a certain threshold.

This formula was tested numerically. It is in very good agreements with the simulations as soon as $S^{\prime}$ is above $\approx 20$. As it decreases, the average of estimated eccentricity increases (which is exactly saying that the bias increases) therefore the approximation of low eccentricities does not hold any more. The value of $S^{\prime}$ can be evaluated keeping in mind that when the linearised model at $e=0$ is poorly conditioned, (matrix M, as defined equations (45) and (46)), then the uncertainty on $k$ and $h$ is higher than given by the simple formula (53) and the $S^{\prime}$ analytical approximation is inoperative.

Appendix D Non Gaussian noise

D.1 Simulations

In this appendix, we show that the non Gaussianity of the noise has a small impact on the quality of the eccentricity estimates.

We generate 1000 realisations of six different types of noises. These have a null mean, are independent, identically distributed and scaled to have a standard deviation $\sigma=1$ (their covariance matrix is the identity). We consider noises that are Gaussian, Student $T$ with 3 and 4 degrees of freedom, uniform, exponential and Poisson. We inject a circular planet on the measurement times of Gl 96 with $K=4\times\sigma$ and uniform $\omega$, $M_{0}$. Ten periods are chosen randomly according to a log normal distribution on the 5 to 500 days interval.

For each type of noise, each of the 10 $\times$ 1000 noise realisations chosen, we compute the error on eccentricity ($|\widehat{e}-e_{t}|$, $e_{t}$ being the true eccentricity and $\widehat{e}$ the estimate), and the root mean square (RMS) of the residuals. The RMS is a s a proxy for the estimated noise level, and thus the width of the error bars. The latter is the maximum likelihood estimate of the noise level for an i.i.d. noise model.

The average values of those on the 10 times 1000 realisations are reported in table 7. Note that we take the square root of the mean squared error (MSE). We find that the mean error on eccentricity for non Gaussian noise is within 2% of the value of the mean error for a Gaussian noise, and the estimated jitter is within 5% of the value of the mean jitter for a Gaussian noise. Only the mean squared error varies by 10% between the Gaussian and Student $T$ distribution with 3 degrees of freedom. Such results are remarkable, since for instance the Poisson noise takes discrete values and is non symmetrical.

The cumulative distribution function of the error and the estimated noise levels (the standard deviation of the residuals) are shown in Fig. 26. Note that the jitter estimates have a slightly greater dispersion for the Student distributions, but as in the case of the average values we do not see striking differences. We simply note that for noise distributions with heavy tails (here Student), as expected, there is a higher fraction of cases where the noise level is severely underestimated or overestimated. We conclude that non Gaussianity does not play a significant role, except for an increased variability of the noise level estimation when the noise distribution has heavy tails.

D.2 Distribution of the bias and the jitter

In this section, we state and prove some mathematical results on the estimate of eccentricity obtained via a linear model. We assume as in section 2.2.2 that the eccentricity estimate is given by the linear model, $\textbf{{y}}=\textbf{{M}}\textbf{{x}}+\mn@boldsymbol{\epsilon}$ where the first four columns of M are such that the $i$-th line is evaluated at observation time $t_{i}$, $\textbf{{M}}_{i}=(\cos(nt_{i})\;\sin(nt_{i})\;\cos(2nt_{i})\;\sin(2nt_{i}))$ and $\textbf{{x}}=(A\;B\;C\;D)^{T}$. The remaining columns M are such that M is of maximal rank. The eccentricity estimate is

[TABLE]

We denote by $\widehat{K}=\sqrt{\widehat{A}^{2}+\widehat{B}^{2}}$. The only assumption on the noise $\mn@boldsymbol{\epsilon}$ is that a vanishing mean and a non degenerate covariance matrix V.

The model just described is supposed to be the true model, which is unknown by the data analyst. We assume that the model with which the analysis is done is $\textbf{{y}}=\textbf{{M}}\textbf{{x}}+\mn@boldsymbol{\epsilon}$ where $\mn@boldsymbol{\epsilon}$ is a Gaussian noise identically distributed of variance $\sigma^{2}$, that is a free parameter. In that model, for a data set $\textbf{{y}}_{0}$ the maximum likelihood estimates of x and $\sigma$ are respectively

[TABLE]

where we denote by I the $N\times N$ identity matrix and by $\textbf{{Q}}=\textbf{{I}}-\textbf{{M}}(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}$.

The question we address is the dependency of the estimate (83) on the noise nature. We show that several quantities relevant to our purposes are only determined by the covariance of the noise. Assertion iv below shows that the bias on $e^{2}$ depends only on the covariance structure of the noise to order 2 in $1/K$. More precisely,

Theorem 1

$\mathbb{E}\{\widehat{\sigma^{2}}\}=\frac{1}{N}\mathrm{tr}(\textbf{{Q}}\textbf{{V}}\textbf{{Q}}^{T})$* where $\mathrm{tr}$ is the sum of the diagonal terms of a matrix (the trace).* 2. 2.

$\mathbb{E}\{\widehat{\textbf{{x}}}\}=\textbf{{x}}_{t}$** 3. 3.

$\mathrm{Cov}\{\widehat{\textbf{{x}}}\}=(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}\textbf{{V}}\textbf{{M}}(\textbf{{M}}^{T}\textbf{{M}})^{-1}$** 4. 4.

$\mathbb{E}\{\widehat{e}^{2}\}-e_{t}^{2}=f(\textbf{{V}})+o\left(\frac{1}{\widehat{K}^{3}}\right)$* .*

Proof

(i) As $\textbf{{y}}_{0}=\textbf{{y}}_{t}+\mn@boldsymbol{\epsilon}$ and $\textbf{{y}}_{t}$ is in the image space of M, we have $\textbf{{Q}}\textbf{{y}}_{0}=\textbf{{Q}}\mn@boldsymbol{\epsilon}$.

The estimate of the noise level is given by eq. (85). Since $\|\textbf{{Q}}\textbf{{y}}_{0}\|^{2}=\textbf{{y}}_{0}^{T}\textbf{{Q}}^{T}\textbf{{Q}}\textbf{{y}}_{0}$ and $\textbf{{Q}}\textbf{{y}}_{0}=\textbf{{Q}}\mn@boldsymbol{\epsilon}$, $\mathbb{E}\{\widehat{\sigma^{2}}\}=\sum_{i,j}(QQ^{T})_{ij}\mathbb{E}\{\epsilon_{i}\epsilon_{j}\}/N=\sum_{i,j}(QQ^{T})_{ij}V_{ij}/N=\frac{1}{N}\mathrm{tr}(\textbf{{Q}}\textbf{{V}}\textbf{{Q}}^{T})$.

(ii) We have \mathbb{E}\{\widehat{\textbf{{x}}}\}=$$\mathbb{E}\{(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}\textbf{{y}}_{0}\}= $(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}(\textbf{{y}}_{t}+\mathbb{E}\{\mn@boldsymbol{\epsilon}\}$ and by hypothesis, $\mathbb{E}\{\epsilon\}=0$, hence the result.

(iii) $\mathrm{Cov}\{\widehat{\textbf{{x}}}\}=$ \mathrm{Cov}\{(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}\textbf{{y}}_{0}\}=$$\mathrm{Cov}\{(\textbf{{M}}^{T}\textbf{{M}})^{-1}\textbf{{M}}^{T}\mn@boldsymbol{\epsilon}\}

(iv) We pose $\widehat{A}=\bar{A}+a$, $\widehat{B}=\bar{B}+b$, $\widehat{C}=\bar{C}+a$, $\widehat{D}=\bar{D}+d$. Denoting by $\bar{K}=\sqrt{\bar{C}+\bar{D}}$. A Taylor expansion of the denominator $\widehat{e^{2}}$ about $\bar{A},\bar{B},\bar{C},\bar{D}$ at order two in $1/\bar{K}$ is

[TABLE]

A simple development shows that $\gamma_{1}=0$ and that $\gamma_{2}=f(\mn@boldsymbol{\mu}^{2})$ where $\mn@boldsymbol{\mu}^{2}$ is the vector of moments of order 2 of x, which is a function of the noise covariance as shown in (iii). For completeness, we give the explicit expression $\gamma_{2}=\sigma_{c}^{2}+\sigma_{d}^{2}+4e_{t}(\cos\psi\cos\phi C_{ac}+\sin\psi\cos\phi C_{cb}+\cos\psi\sin\phi C_{ad}+\sin\psi\sin\phi C_{db})+e_{t}^{2}(\sigma_{a}^{2}+\sigma_{b}^{2})$ where $C_{ij}$ is the covariance of $\widehat{I}$ and $\widehat{J}$ and $\bar{A}=\bar{K}\cos\phi$, $\bar{B}=\bar{K}\sin\phi$ $\bar{C}=\bar{e}\bar{K}\cos\phi$, $\bar{C}=\bar{e}\bar{K}\sin\phi$.

Appendix E Residual analysis

In this section we compute the law followed by the residuals of a linear least-square fit. Let us suppose that we have a model

[TABLE]

where y is a vector of $N$ observations, modelled as a linear combination of the column of the $N\times p$ matrix $A$, and $\mn@boldsymbol{\epsilon}$ is a Gaussian noise of covariance matrix $\textbf{{V}}=:\textbf{{W}}^{-1}$. Assuming V and A are known, the least square estimate of y is $\widehat{\textbf{{y}}}=\textbf{{A}}(\textbf{{A}}^{T}\textbf{{W}}\textbf{{A}})^{-1}\textbf{{A}}^{T}\textbf{{W}}\textbf{{y}}$ Therefore

[TABLE]

where $\textbf{{I}}_{N}$ is the identity matrix of size $N$ and $\textbf{{B}}:=\textbf{{A }}(\textbf{{A}}^{T}\textbf{{W}}\textbf{{A}})^{-1}\textbf{{A}}^{T}\textbf{{W}}$. The quantity $\mathit{r}_{W}$, being a product of a matrix ($\textbf{{W}}^{1/2}(\textbf{{I}}_{n}-\textbf{{B}})$) with a Gaussian random variable of covariance V has a covariance U

[TABLE]

where $\textbf{{C}}=\textbf{{W}}^{1/2}\textbf{{A}}$. This notation is convenient because it shows clearly that $\textbf{{P}}=\textbf{{C}}(\textbf{{C}}^{T}\textbf{{C}})^{-1}\textbf{{C}}^{T}$ is a projection matrix on the space generated by the columns of C. Finally

[TABLE]

Is a projection on the space orthogonal to the one generated by C columns. Therefore, there exists an orthonormal matrix Q such that

[TABLE]

where $\textbf{{J}}_{p}$ is a diagonal matrix whose first $p$ elements are zero and the others are equal to one. Finally, let us remark that $\textbf{{r}}_{QW}:=\textbf{{Q}}^{T}\textbf{{r}}_{W}$ has a covariance matrix $\textbf{{Q}}^{T}\textbf{{UQ}}=\textbf{{J}}_{p}$, which shows the claim of section 4, $\textbf{{Q}}^{T}\textbf{{W}}^{1/2}(\textbf{{y}}-\widehat{\textbf{{y}}})$ has $p$ 0 components and the others are Gaussian variables of mean [math] and variance 1.

Let us finally note that the covariance matrix $\textbf{{U}}=\textbf{{I}}_{N}-\textbf{{P}}$ of $\textbf{{r}}_{W}=\textbf{{W}}^{1/2}(\textbf{{y}}-\widehat{\textbf{{y}}})$ will be close to identity if there are many more observations than parameters. This explains why the weighted residuals $\textbf{{r}}_{W}$ almost behaves like independent Gaussian variables and, for instance, why plotting $\textbf{{r}}_{W}(t_{i})-\textbf{{r}}_{W}(t_{j})$ as a function of $t_{i}-t_{j}$ gives hints on the correlations.