Evidence for at least three planet candidates orbiting HD20794

TL;DR

This study uses advanced noise modeling on HARPS data to identify at least three potential Earth-like planets orbiting HD20794, including a newly reported super-Earth at 147 days and a candidate in the habitable zone.

Contribution

The paper introduces a novel noise diagnosis method using differential radial velocities, enabling improved detection of low-mass planet candidates around Sun-like stars.

Findings

Detected four signals at 18, 89, 147, and 330 days.

First report of a super-Earth at 147 days.

Evidence for a potential habitable zone planet at 330 days.

Abstract

We explore the feasibility of detecting Earth analogs around Sun-like stars using the radial velocity method by investigating one of the largest radial velocities datasets for the one of the most stable radial-velocity stars HD20794. We proceed by disentangling the Keplerian signals from correlated noise and activity-induced variability. We diagnose the noise using the differences between radial velocities measured at different wavelength ranges, so-called "differential radial velocities". We apply this method to the radial velocities measured by HARPS, and identify four signals at 18, 89, 147 and 330 d. The two signals at periods of 18 and 89 d are previously reported and are better quantified in this work. The signal at a period of about 147 d is reported for the first time, and corresponds to a super-Earth with a minimum mass of 4.59 Earth mass located 0.51 AU from HD20794. We also find a significant signal at a period of about 330 d corresponding to a super-Earth or Neptune in the habitable zone. Since this signal is close to the annual sampling period and significant periodogram power in some noise proxies are found close to this signal, further observations and analyses are required to confirm it. The analyses of the eccentricity and consistency of signals provide weak evidence for the existence of the previously reported 43 d signal and a new signal at a period of about 11.9 d with a semi amplitude of 0.4 m/s. We find that the detection of a number of signals with radial velocity variations around 0.5\,m/s likely caused by low mass planet candidates demonstrates the important role of noise modeling in searching for Earth analogs.

Figures

Tables

Table 1. Table 1: The prior distributions of model parameters. The unit of c k subscript 𝑐 𝑘 c_{k} and d m subscript 𝑑 𝑚 d_{m} is m/s. The maximum and minimum time of the data are denoted by t max subscript 𝑡 max t_{\textrm{max}} and t min subscript 𝑡 min t_{\textrm{min}} , respectively. The maximum RV with respect to the mean is denoted as | v − v ¯ | max subscript 𝑣 ¯ 𝑣 max |v-\bar{v}|_{\rm{max}} . The parameter characterizing the linear dependence of RV on activity indices is c k subscript 𝑐 𝑘 c_{k} , where k 𝑘 k denotes various indices.

Parameter	Unit	Prior distribution	Minimum	Maximum
Each Keplerian signal
$K_j$	m/s	$1/\left(K_{\text{max}}-K_{\text{min}}\right)$	0	$2{\left|v-\overline{v}\right|}_{\text{max}}$
$P_j$	day	$P_j^{-1}/\log \left(P_{\text{max}}/P_{\text{min}}\right)$	1	$t_{\text{max}}-t_{\text{min}}$
$e_j$	—	$𝒩(0,0.1)$	0	1
${\omega }_j$	rad	$1/\left(2\pi \right)$	0	$2\pi$
$M_{0j}$	rad	$1/\left(2\pi \right)$	0	$2\pi$
Linear trend and jitter
$a$	m s-1yr-1	$1/\left(a_{\text{max}}-a_{\text{min}}\right)$	$-365.24K_{\text{max}}/P_{\text{max}}$	$365.24K_{\text{max}}/P_{\text{max}}$
$b$	m/s	$1/\left(b_{\text{max}}-b_{\text{min}}\right)$	$-K_{\text{max}}$	$K_{\text{max}}$
$s_J$	m/s	$1/\left(s_{J\max }-s_{J\min }\right)$	0	$K_{\max }$
Moving average
$w$	—	$1/\left(w_{\text{max}}-w_{\text{min}}\right)$	-1	1
$\tau$	day	$P_j^{-1}/\log \left({\tau }_{\text{max}}/{\tau }_{\min }\right)$	$1/\left(t_{\max }-t_{\min }\right)$	1
Activity indexes and differential RVs
$c_k$	m/s	$1/\left(c_{k\text{max}}-c_{k\text{min}}\right)$	$-c_{k\text{max}}$	$K_{\text{max}}/\left(I_{X\text{max}}-I_{X\text{min}}\right)$
$d_m$ ($m\in \{1,\dots ,N_{\lambda }-1\}$)	m/s	$1/\left(d_{m\max }-d_{m\text{min}}\right)$	$-d_{m\max }$	$K_{\text{max}}/\left(D_{m\max }-D_{m\min }\right)$

Table 2. Table 2: The logarithmic BIC-estimated BFs for noise models with respect to the model without any MA components or differential RVs. For all models, the BIS, FWHM and S-index are linearly combined with differential RVs.

$N_D$	MA(0)	MA(1)	MA(2)	MA(3)	MA(4)	MA(5)	MA(6)
0	0	1769	1969	2006	2011	2008	2008
2	151	1963	2219	2288	2299	2299	2299
5	238	1985	2241	2312	2325	2324	2323
8	236	1976	2233	2305	2318	2317	2316

Table 3. Table 3: The logarithmic BIC-estimated BFs for models with various numbers of planets for Keplerian and circular solutions.

Number of planets	1	2	3	4	5	6
Period in Keplerian solution	332.69	18.33	88.98	147.08	11.86	43.14
$\log \left(\mathrm{BF}\right)$	38	72	96	134	137	135
Period in circular solution	337.26	18.32	86.80	147.33	43.33	11.29
$\log \left(\mathrm{BF}\right)$	46	83	95	121	145	147

Table 4. Table 4: The maximum a posteriori (MAP) estimation of the parameters for the noise model (lower rows) and the five signals (upper rows) detected in the TERRA-reduced HARPS data for HD20794. The parameters are estimated in the 6-planet Keplerian solution. The uncertainties of parameters are represented by the values determined at 1% and 99% of the cumulative posterior density. We estimate the minimum planetary mass and semi-major axis using a stellar mass of 0.813  M ⊙ − 0.012 + 0.018 superscript subscript subscript 𝑀 direct-product 0.012 0.018 {}^{+0.018}_{-0.012}M_{\odot} (Ramírez et al. 2013 ) . Among the parameters of the noise model, c 1 subscript 𝑐 1 c_{1} , c 2 subscript 𝑐 2 c_{2} and c 3 subscript 𝑐 3 c_{3} are linear coefficients for S-index, BIS and FWHM, respectively.

Parameters	HD20794 b	HD20794 d	HD20794 e	HD20794 f?	HD20794 g?	HD20794 c?
$P$ (d)	18.33 [18.31, 18.34]	88.90 [88.49, 89.27]	147.02 [146.11, 148.45]	331.41 [328.40, 336.49]	11.86 [11.84, 11.87]	43.17 [43.07, 43.29]
$K$ (m/s)	0.81 [0.57, 0.81]	0.60 [0.42, 0.70]	0.69 [0.56, 0.84]	1.09 [0.94, 1.28]	0.34 [0.24, 0.50]	0.53 [0.36, 0.64]
$e$	0.27 [0.05, 0.31]	0.25 [0.04, 0.41]	0.29 [0.11, 0.43]	0.05 [0.00, 0.11]	0.20 [0.01, 0.35]	0.17 [0.01, 0.27]
$\omega$ (rad)	6.74 [5.75, 7.22]	4.41 [3.64, 5.10]	4.68 [4.11, 5.52]	3.24 [0.91, 5.64]	0.77 [-1.47, 2.05]	0.55 [-1.60, 1.91]
$M_0$ (rad)	2.26 [1.56, 3.07]	4.69 [3.89, 5.60]	2.52 [1.76, 3.46]	4.02 [1.00, 5.72]	2.74 [0.96, 4.74]	-1.67 [-2.87, 1.68]
$m\sin i$ ($M_{\oplus }$)	2.82 [2.02, 2.92]	3.52 [2.51, 4.10]	4.77 [3.91, 5.73]	10.26 [8.79, 12.15]	1.03 [0.73, 1.52]	2.52 [1.69, 3.04]
$a$ (au)	0.127 [0.126, 0.128]	0.364 [0.360, 0.368]	0.509 [0.503, 0.515]	0.875 [0.865, 0.886]	0.095 [0.094, 0.096]	0.225 [0.222, 0.227]
$a$ (m s-1 yr-1)	$b$ (m s-1)	$c_1$ (m s-1)	$c_2$ (m s-1)	$c_3$ (m s-1)	$s_J$ (m s-1)	$w_1$
-0.52 [-0.64, -0.47]	-1.69 [-2.02, -1.12]	-0.07 [-0.13, 0.01]	0.02 [-0.01, 0.06]	0.93 [0.82, 1.06]	0.89 [0.87, 0.92]	0.41 [0.37, 0.44]
$w_2$	$w_3$	$w_4$	$\ln \left(\tau \right)$ (day)	$d_1$ (m s-1)	$d_2$ (m s-1)	$d_3$ (m s-1)
0.23 [0.20, 0.28]	0.19 [0.14, 0.23]	0.14 [0.07, 0.16]	-2.06 [-2.37, -1.75]	-0.15 [-0.20, -0.13]	-0.35 [-0.40, -0.31]	-0.40 [-0.46, -0.36]
$d_4$ (m s-1)	$d_5$ (m s-1)
-0.36 [-0.42, -0.31]	-0.23 [-0.28, -0.19]

Equations

$$D(t_i,{\lambda }_j)\equiv v(t_i,{\lambda }_{j+1})-v(t_i,{\lambda }_j)=\Psi (t_i,{\lambda }_{j+1})-\Psi (t_i,{\lambda }_j)\text{leavevmode},$$

$$\hat{\Psi }(t_i,{\lambda }_j)=a({\lambda }_j)t_i+b({\lambda }_j)+\sum _k{c}_k({\lambda }_j)I_k+\sum _{m=1}^{N_{\lambda }-1}{d}_{m}D(t_i,{\lambda }_m)\text{leavevmode},$$

$${\hat{v}}_b(t_i,{\lambda }_j)$$

$$f_k(t_i)$$

$$\hat{v}(t_i,{\lambda }_j)={\hat{v}}_b(t_i,{\lambda }_j)+\sum _k{w}_k({\lambda }_j)\exp [-|t_i-t_{i-k}|/\tau ({\lambda }_j)]{ϵ}_{i-k},$$

$${\mathcal{L}}_j\equiv P(v_j|\mathbf{\theta },M)=\prod _i\frac{1}{\sqrt{2\pi [{\sigma }_{i,j}^2+s_J({\lambda }_j{)}^2]}}\exp \left[-\frac{(\hat{v}(t_i,{\lambda }_j)-v_{i,j}{)}^2}{2({\sigma }_{i,j}^2+s_J({\lambda }_j{)}^2)}\right]\text{leavevmode},$$

Full text

*Methods. * We proceed by disentangling the Keplerian signals from correlated noise and activity-induced variability. We diagnose the noise using the differences between radial velocities measured at different wavelength ranges, so-called “differential radial velocities”.

*Results. *We apply this method to the radial velocities measured by HARPS, and identify four signals at 18, 89, 147 and 330 d. The two signals at periods of 18 and 89 d are previously reported and are better quantified in this work. The signal at a period of about 147 d is reported for the first time, and corresponds to a super-Earth with a minimum mass of 4.59 Earth mass located 0.51 AU from HD20794. We also find a significant signal at a period of about 330 d corresponding to a super-Earth or Neptune in the habitable zone. Since this signal is close to the annual sampling period and significant periodogram power in some noise proxies are found close to this signal, further observations and analyses are required to confirm it. The analyses of the eccentricity and consistency of signals provide weak evidence for the existence of the previously reported 43 d signal and a new signal at a period of about 11.9 d with a semi amplitude of 0.4 m/s.

*Conclusions. *We find that the detection of a number of signals with radial velocity variations around 0.5 m/s likely caused by low mass planet candidates demonstrates the important role of noise modeling in searching for Earth analogs.

††offprints: Fabo Feng,

1 Introduction

High precision spectrometers enable us to find super-Earths by measuring the Doppler shift of the stellar spectra caused by the periodic perturbations of exoplanets orbiting the target stars. Unlike the transit method which relies on a particular orbital orientation, the radial velocity (RV) method can provide useful information about planets around all stars. The recent detection of the nearest “habitable-zone” planet candidate Proxima Centauri b demonstrates the important role of the RV technique in detecting Earth analogs (Anglada-Escudé et al. 2016).

However, the stellar spectrum is contaminated by various noise sources such as stellar activity, instrumental noise and non-perfect data reduction. An analysis of high-cadence spectroscopy of M dwarfs shows that the noise floor (about 1 m/s) of current Doppler surveys has to a large extent an instrumental origin (Berdiñas et al. 2016). Our analysis is also concerned with demonstration of the strong dependence of activity-induced noise on wavelength. The minimization of both instrumental and stellar noise requires a noise model accounting for noise correlated in time and wavelength.

Considering that complex models typically lead to false negatives while simple models lead to false positives (Feng et al. 2016), we follow the Goldilocks principle introduced in Feng et al. (2016) to select the best noise model for a given RV data set. In other words, we devise a noise model framework to diagnose the noise in a given data set. In this work, we apply this noise model framework to the HARPS RV measurements of HD20794.

HD20794 is a high-velocity and metal-deficient G8 star with a mass of $0.813^{+0.018}_{-0.012}$ $M_{\odot}$ (Ramírez et al. 2013). It has been reported to host at least 3 planets and a dust disk (Pepe et al. 2011; Kennedy et al. 2015). The postulated planets are likely super-Earths with orbital periods of 18.1, 40.1 and 90.3 d (Pepe et al. 2011) (hereafter P11). These results of P11 were obtained using periodogram analysis which assumes white RV noise and circular planetary orbits. Although HD20794 is known to have planets, it has long been chosen as a special RV reference target due to its notable stability over long time baselines found by the Southern radial velocity programmes of the Anglo-Australian Planet Search (e.g., Fig. 2 of Butler et al. 2001) and its continuous presence on the HARPS Guaranteed Time Observation list111e.g. http://www.eso.org/sci/observing/teles-alloc/gto.html.

This paper is structured as follows. We introduce the RV data in section 2. In section 3, we introduce differential RVs to model wavelength-dependent noise and choose the so-called Goldilocks noise model in the Bayesian framework. In section 4, we apply our models to the data and report the results. Finally, we discuss and conclude in section 5.

2 HARPS Doppler measurements for HD20794

In the European Southern Observatory archive, there are 5150 publicly available RVs measured by HARPS from September 2003 to September 2013 for HD20794 from programmes Mayor, 60.A-9036, 072.C-0488, 083.C-1001, 084.C-0229, 086.C-0230, 087.C-0990, 088.C-0011, 089.C-0050, 090.C-0849, 091.C-0936. We process the Doppler measurements using the TERRA algorithm (Anglada-Escudé & Butler 2012), which generates RVs for 72 echelle orders. The stellar activity is recorded by the following indices: the S-index measured from the Ca II H&K emission, line bisector shapes (BIS), and the width of the spectral lines (FWHM). We exclude the points with RVs or activity indices which deviate from the mean by more than 5-$\sigma$. We do this iteratively to omit any visible outliers. Moreover, we also remove the data obtained before JD2453500 due to the large dispersion in the activity indices, FWHM in particular. Such a feature is also observed in the $\tau$ Ceti data (Feng et al. 2017b), although the cause for this may be different. Finally, we obtain a data set of 4882 RVs together with activity indices. These data and the discarded epochs are shown in Fig. 1. We see a great variation of FWHM before JD2453500, and this feature is also found in the data measured for $\tau$ Ceti (Feng et al. 2017b). The data is available in http://star-www.herts.ac.uk/~ffeng/HD20794_supplementary/data.

We use the RV differences between the 72 echelle orders to remove wavelength-dependent noise. We divide the 72 spectral orders into groups, and average the RVs in each group to generate the so-called “aperture data sets”. For $n$ groups of spectral orders, we generate $n$AP$j$ aperture data sets, where $j$ is a natural number not larger than $n$. We define the RV differences between aperture data sets as differential RVs. We denote them by $n$AP$j-i$, where $n$ is the number of groups of spectral orders, and $j$ and $i$ denote different aperture data sets. For $n$ independent aperture data sets, there are $n-1$ independent differential RVs in total. The schematic of this data reduction process is shown in Fig. 2.

3 Noise model framework and Bayesian inference

3.1 Wavelength-dependent noise models

The RV variation induced by stellar activity and rotation is typically wavelength-dependent (Tuomi et al. 2013a). Thus the average of RVs measured at different wavelength ranges (or spectral orders) is biased due to the ignorance of such a dependence. To remove such bias, we introduce “differential RVs” to weight the spectral orders a posteriori. Differential RVs are RV differences between different spectral orders (or wavelength ranges) and thus do not contain Keplerian variation since Keplerian variation is color-independent. Following the Goldilocks principle of noise modeling introduced by Feng et al. (2016), we reduce the number of different RVs by deriving them from averaged spectral orders, so-called aperture data sets (see Fig. 2).

Then the differential RVs are

[TABLE]

where $v(t_{i},\lambda_{j})$ is the RV measured at time $t_{i}$ within a wavelength range centered at $\lambda_{j}$. We define the noise component of the RV model as

[TABLE]

where $N_{\lambda}$ is the number of wavelength ranges or aperture divisions and thus $N_{D}\equiv N_{\lambda}-1$ is the number of independent differential RVs, $d_{m}$ characterizes the linear dependence of RV noise on the $m^{\rm th}$ differential RV, $a$ is the linear acceleration caused by stellar companions or activity cycles and $b$ is the reference velocity. The linear dependence of RVs on activity index $I_{k}$ is parameterized by constant $c_{k}$. We consider that $a$ and $b$ can be wavelength dependent, because besides fitting for trends due to companions, those parameters can also absorbs any drift due to long term activity not correlated with activity indices or long term instrumental noise variation. We linearly combine $\hat{\Psi}(t_{i},\lambda_{j})$ with $n$ Keplerian components to define the basic RV model, which is

[TABLE]

where $f_{k}(t_{i})$ is the $k^{\mathrm{th}}$ Keplerian component, and $\hat{\Psi}$ is the noise component. In the Keplerian function, $K_{k}$, $\omega_{k}$, $\nu_{k}$, $e_{k}$ are the amplitude, longitude of periastron, true anomaly and eccentricity for the $k^{\mathrm{th}}$ Keplerian signal.

Since the moving average (MA) model can efficient account for red noise in RVs and avoid false negatives (Feng et al. 2016), we model the red noise using a general MA model with exponential smoothing which is

[TABLE]

where $w_{k}$ and $\tau$ are the amplitude and time scale of the moving average, and $\epsilon_{i-k}$ is $v(t_{i-k},\lambda_{j})-\hat{v}_{b}(t_{i-k},\lambda_{j})$. Hereafter, we define “nP+MA(q)+mD” as the $n$-planet model combined with $q^{\rm{th}}$-order moving average and $m$ differential RVs deriving from $m+1$ aperture data sets (see Fig. 2). The white noise model is $\hat{\Psi}(t_{i},\lambda_{j})$, denoted by MA(0).

We also account for a white excess noise with an amplitude of $s_{J}$ in the likelihood which is

[TABLE]

where $\sigma_{i,j}$ is the measurement noise at time $t_{i}$ in the jth aperture data set with wavelength range centered at $\lambda_{j}$, $s_{J}(\lambda_{j})$ is the jitter level, and $v_{i,j}$ is the RV measured at time $t_{i}$ in the jth aperture data set. Following Feng et al. (2016), we adopt a Gaussian prior distribution centered at zero and with a standard deviation of 0.1 for eccentricity, logarithmic uniform distributions for period $P$ and correlation time scale $\tau$, and uniform distributions for other parameters. The specific prior distributions are given in Table 1.

In almost all previous research, only the averaged data or 1AP1 is analyzed without accounting for wavelength-dependent noise. We justify the usage of differential RVs in the removal of wavelength-dependent noise by applying the model defined in Eqn. 4 and 5 to the 3AP1, 3AP2, 3AP3 and 1AP1 data sets. We model these data sets using the noise models of 0P+MA(1)+0D, 0P+MA(1)+2D, 0P+MA(0)+0D, 0P+MA(0)+2D. We calculate the generalized periodogram with floating trend (GLST; Feng et al. 2017a) for each data set and each model to visualize the consistency of signals in wavelength. Following Cumming et al. (1999), the GLST is normalized with respect to the residual variance and the FAP is calculated accordingly.

The GLSTs for the 3AP1, 3AP2, 3AP3 and 1AP1 data sets and their residuals after subtracting the best-fit noise model are shown in Fig. 3. We observe significant differences between GLSTs for 3AP1, 3AP2, 3AP3 and 1AP1, indicating significant wavelength-dependent noise in the RV data. On the contrary, the residuals of all RV data sets have similar GLSTs after subtracting the best-fit 0P+MA(0)+2D or 0P+MA(1)+2D from the data. Such consistency is not sensitive to the change of time-correlated noise component (i.e. MA(0) and MA(1)), as seen in the rightmost two columns. Therefore differential RVs are essential in removing wavelength-dependent noise.

In the following sections, we will focus our analysis on the 1AP1 data set because it has a higher signal to noise ratio with respect to the other aperture data sets. As seen in Fig. 3, the wavelength-dependent noise in 1AP1 is properly removed by differential RVs a posteriori.

3.2 Goldilocks noise model

We apply the adaptive Metropolis Markov Chain Monte Carlo (AM) introduced by Haario et al. (2001) to sample the posterior. Specifically, we use tempered (or hot) AM chains to identify potential signals, and then use un-tempered (or cold) chains to constrain signals. Similar methods have been introduced by Gregory (2011). The hot chain is defined as a chain moving according to the modified posterior $P({\mathcal{D}}|\theta,M)^{\beta}P(\theta|M)$ with $0<\beta\leq 1$, where $P({\mathcal{D}}|\theta,M)$ and $P(\theta|M)$ are the likelihood and prior of parameters $\theta$ for model $M$ and data ${\mathcal{D}}$. We use four parallel chains for inference if the posterior samples drawn by these chains converge to a stationary distribution. In other words, the ratio of variances between and within these chains should be less than 1.1, which is the so-called “Gelman-Rubin criteria” (Gelman & Rubin 1992). We typically use the AM algorithm to draw a few million posterior samples for model and/or parameter inference. The detailed results are available online: http://star-www.herts.ac.uk/~ffeng/HD20794_supplementary/results.

From the likelihoods of the posterior samples drawn by cold chains, we estimate the Bayes factor (BF) using the Bayesian information criterion (BIC). Following Feng et al. (2016), we regard model A as favored over model B if the logarithmic BF of the former with respect to the latter is larger than 5. To confirm a signal, we also adopt the criteria used by Tuomi (2012) that the period should be constrained from above and below in the posterior distribution of the period. In other words, $P(P_{k}|{\mathcal{D}},\leavevmode\nobreak\ M)$ converges to a stationary distribution.

Following the Goldilocks principle introduced by Feng et al. (2016), we apply the BF criterion to determine which noise model to use for removing noise. We generate a set of noise models by setting $q\in\{0,1,2,3,4,5,6\}$ and $N_{D}\in\{0,2,5,8\}$, and apply them to the data. The results of the model comparison are shown in Table 2. We find that the fourth order MA combined with five differential RVs is the best choice. Apart from $q$ and $N_{D}$, all other model parameters are free. Specifically there are 2 parameters for the trend, 1 for jitter, 3 for activity indices, 5 for differential RVs and 5 for the MA model. There are 16 free parameters in total in the 0P+MA(4)+5D model. In the following section, we analyze the HARPS data using this model.

4 Keplerian candidates

We apply both Keplerian and sinusoidal functions to model the RV variations caused by planets, which are called “Keplerian solution” and “circular solution”, respectively. We identify six signals in both solutions, and report the BFs for them in Table 3. Although the $\sim$43 d signal identified in the Keplerian solution does not pass the BF threshold of 150, it satisfies the criteria for the circular solution. Moreover, without including the 43 d signal into the model, the eccentricity of the 89 d signal would be as high as 0.4. The 43 d and 89 d signals are approximately in 1:2 resonance, which is not rare for exoplanets (Steffen & Hwang 2015). Although the signal at a period of 11 d does not pass the BF threshold in both solutions, it is strong in the Keplerian solution, and increases the BF by at least one order of magnitude. In addition, it is visible in the periodograms for all observation seasons, as we will see later. The periods in the circular solution are different from the ones in the Keplerian solution because of the assumption of periods in the circular solution.

To test whether the identified signals are caused by stellar and instrumental noise, we show the GLSTs for the data, the activity indices and differential RVs in Fig. 4. In the figure, we find that the RVs are strongly correlated with the FWHM and differential RVs, demonstrating the essential role of these noise proxies in removing correlated noise. We also find that periods around the $\sim$330 d signal are strong in noise proxies such as FWHM and 6AP5-4, though these periods typically deviate from 330 d by more than 20 d. Such approximate overlaps are expected since the 330 d signal is close to the annual sampling frequency. Moreover, the probability of a random overlap between a signal and significant powers in the GLSTs of noise proxies increase with the number of proxies. Considering the minimization of correlated noise using noise proxies in the model, the 330 d signal is too strong to be caused by pure noise.

We also calculate the GLSTs for the data subtracted by the optimal prediction of models with different numbers of Keplerian components, and show them in Fig. 5. We see that all identified signals correspond to certain peaks in the periodogram for the original RVs. The subtraction of the noise components makes some signals more significant. We find that the signal identified by posterior sampling is not always the most significant one shown in the residual periodogram. Moreover, the 43 and 11 d signals are not strong in the corresponding residual periodograms, indicating the limitation of residual periodograms in identifying Keplerian signals. The residuals after subtracting the predictions of n-planet model are shown in Fig. 8 in the appendix. The correlation between the RV residuals and noise proxies are shown in Fig. 9.

To show the consistency of the detected signals, we make a Bayesian periodogram for data within a moving time window. To make this periodogram, we model the RVs using a combination of sinusoidal functions and a linear trend, and analytically marginalize the likelihood over the amplitudes of sinusoidal functions and the trend parameters which are $a$ and $b$ in Eqn. 2. Hence this new periodogram is called “marginalized likelihood periodogram” (MLP; Feng et al. 2017a). In addition to trend parameter $a$, the MLP also optimize $b$, and thus extend the Bayesian generalized Lomb-Scargle periodogram (Mortier et al. 2015)222The relevant code for MLP and GLST is available in Github: https://github.com/phillippro/agatha, based on which a web app is also developed and linked at http://www.agatha.herts.ac.uk. To make a 2D periodogram, we move the time window in 100 steps to cover the whole data set, and compute the MLP for each step to construct a moving periodogram. We apply this method to the RVs from which the noise component in the 6-planet model prediction for the Keplerian solution has been subtracted. The moving periodograms with 1000 d and 2000 d time windows together with the data are shown in Fig. 6. Since more data is needed to cover the phase of long period signals, we recommend the 2000 d time window for testing the consistency of the 330 d and 147 d signals. We calculate the MLP for each time step, and scale the logarithmic marginalized likelihood (ML) to be ${\rm RML}\equiv({\rm ML}-\overline{\rm ML})/({\rm ML}_{\rm max}-\overline{\rm ML})$ where $\overline{\rm ML}$ and ${\rm ML}_{\rm max}$ are the mean and maximum ML. The RMLs are encoded by colors to show the consistency of signals in time rather than to estimate the significance of signals.

As we can see, the 11 d signal together with its aliases, though weak, are visible over the whole time span. The 147 and 330 d signals are consistently identified in the 2000 d moving periodogram. Notably the significance of these signals may vary with time because long period signals are sensitive to data sampling and size, which typically vary with time. There are also strong powers around the harmonics of the 330 d signal at a period of about 680 d and the annual alias of the 147 d signal at a period of about 250 d. Considering that the 330 d signal is also strong in some noise proxies, further observations and analyses are required to explore the nature of this long period signal.

The 147 d signal is strong in many observation seasons despite the slight variation of its period. The high periodogram power around 110 d is probably a combination of the aliases of 89 and 147 d since the subtraction of the 89 d signal would reduce its significance. The signals around 89 d, 43 d and 18 d are also consistently identified over the whole time span. These signals are previously found by P11 in the circular solutions. Notably the periods of $\sim$43 and 147 d slightly deviate from the optimized values shown in the periodogram because the periodogram assumes zero eccentricity. In addition to the identified signals, the periodogram maps in Fig. 6 also show annual alias of 330 d around 174 d, and annual aliases of 89 d around 113 and 70 d. The annual aliases of 43, 18 and 11 d are also visible. Therefore we interpret all strong signals based on the moving periodograms which visualize the genuine signals as well as their aliases.

Based on the above analyses, we interpret the signals around 18, 89, 147 and 330 d as planet candidates, though investigations of the 330 d signal is necessary to fully confirm its Keplerian nature. Since there is only weak evidence for the signals at periods of 11 and 43 d, more data and further analyses are required to confirm these potential candidates. We show the phase curves for these six signals in Fig. 7. We see that the 6-planet model well fits the data, supporting the existence of 6 signals despite considerable eccentricity for signals around 18, 43 and 147 d.

In Table 4, we report the parameters of the noise model and the four planet candidates, which are estimated in the 6-planet Keplerian solution. In the table, we see high MAP values of linear parameters such as $c_{3}$ and $d_{3}$333These parameters measure the contribution of noise RV variations and are in units of m/s because the proxies are normalized. On the other hand, the Pearson coefficients shown in Fig. 9 measure the correlation between RVs and proxies., suggesting the existence of correlated noise as high as 1 m/s. This correlated noise would prevent any reliable detections of sub-meter signals if not minimized through the application of noise proxies such as activity indices and differential RVs. However, the application of noise proxies may introduce extra noise in the fitting (Feng et al. 2017b) due to an oversimplification of the complex/nonlinear dependence of RVs on proxies by linear functions. This is also part of the reason why we do not interpret the weak signals at periods of 11 and 43 d as planet candidates. On the other hand, the MA model plays an important role in removing time-correlated noise, which is probably intra-night RV variability caused by instrumental or reduction-process effect (Berdiñas et al. 2016).

In Table 4, the signal at a period of 147 d corresponds to a super-Earth perturbing HD20794 with a mean semi-amplitude444The mean semi-amplitude is the semi-amplitude $K$ multiplied by $1-e$. of about 0.5 m/s. The candidate with 89 d period has a minimum mass lower than the candidate around 90 d reported by P11 probably because we have accounted for correlated noise which might have been interpreted as signals by P11. The 11 d signal caused a radial velocity variation of about 0.35 m/s. Regardless of whether this signal is Keplerian or not, the detection of such a weak signal demonstrates the important role of noise modeling in detecting low mass planets using the radial velocity technique.

According to the results of RV challenge, the MA model in combination of Bayesian methods can identify signals with a $K/N$ ratio as low as 5 without announcing a false positive (Dumusque et al. 2017). The K/N ratio for signal with semi-amplitude of $K$ is defined as $K/N\equiv K/RV_{\rm rms}\times\sqrt{N_{\rm obs}}$, where $RV_{\rm rms}$ is the standard deviation of RVs after removing the best-fit trend and correlation with noise proxies, and $N_{\rm obs}$ is the number of observations. If we consider measurements in a 15 min bin as an independent observation (Mayor et al. 2003; O’Toole et al. 2008)555The RVs were measured with high cadence and are thus highly correlated in time due to stellar oscillations and systematic errors (Teixeira et al. 2009). It is necessary to bin the data to calculate the number of independent observations. , we get 713 independent observations and $RV_{\rm rms}=1.61$ m/s. Then the $K/N$ ratios are 12, 10, 10, 19, 6.7 and 9.6 for signals at periods of 18, 89, 147, 330, 11 and 43 d. Hence the sub-meter signals we have identified for HD20794 b, d, and e are well above the $K/N=7.5$ threshold of reliable detections of RV signals. Even the smallest signal at 11 d is close to this threshold and is higher than the $K/N=5$ limit reached by our team.

However, the nominal eccentricities of HD20794 b, d, and e are different from zero, probably leading to dynamical instability of the system according to the Lagrange stability analyses (Barnes & Greenberg 2006). But we find that the system is stable if lower eccentricity within the uncertainty interval is adopted for HD20794 b, d and e. The high eccentricity reported in Table 4 might arise from instrumental noise (Feng et al. 2017b). Specifically, the instrumental noise may cause short-term RV variations which favor high eccentricity solutions. The high eccentricity may also be caused by the intrinsic bias in estimating non-negative eccentricity (Zakamska et al. 2011). Therefore the proposed candidates probably have orbits more circular than those reported in this work.

Only the candidate with a period of about 330 d is located in the habitable zone (Kopparapu et al. 2013), although the confirmation of this candidate requires further observations and analyses. Its habitability would be influenced by the frequent impacts of objects from the massive circumstellar disc (Kennedy et al. 2015). All the other candidates are located too close to the host star to allow the existence of liquid water, and thus cannot be considered habitable.

5 Discussions and conclusions

We analyze the HARPS data of HD20794 in the Bayesian framework. We find strong dependence of the RV noise on wavelengths. This wavelength-dependent noise cannot be removed by averaging all spectral orders as previous studies did. To deal with this noise, we use differential RVs to weight the spectral orders a posteriori. We apply this method to data sets measured within different wavelength ranges, and find that differential RVs efficiently remove wavelength-dependent noise. Therefore we propose a combination of the MA model and differential RVs to remove the time and wavelength-dependent noise in RV datasets.

By modeling the RV noise correlated in time and wavelength for HD20794, we identify three firm Keplerian candidates at periods of about 18, 89 and 147 d. The signal at a period of 330 d is probably Keplerian, although similar periods have strong powers in some proxies. There is also weak evidence for the existence of signals at periods of 43 and 11 d. While the 43 d signal does not pass the Bayes factor criterion, the inclusion of it in the model can reduce the eccentricity of the 89 d signal. Although the 11 d signal only increases the Bayes factor by a factor of about 20, it consistently appears in different data chunks.

Thus by applying noise modeling to the P11 and newer HARPS data, we confirm the two candidates at periods of 18 and 90 d reported by P11 and quantify them better. But we are suspicious about the existence of the reported 40 d candidate because it does not pass the signal detection criteria. We also find a new planet candidates, HD20794 e, at a period of about 147 d, corresponding to a super-Earth. The signal at a period of about 330 d probably correspond to a Neptune located in the habitable zone of HD20794. The estimated orbital eccentricity of HD20794 b, d and e is larger than 0.2.

This considerable eccentricity is probably caused by instrumental noise and the fitting bias, as concluded by (Feng et al. 2017b) and Zakamska et al. (2011). Notably such eccentric solutions are also found in the HARPS data of $\tau$ Ceti (Feng et al. 2017b). This indicates a connection between high eccentricity and significant red noise in large RV datasets of bright stars that are observed with high cadence. Thus it is a viable hypothesis that the four planet candidates have eccentricities low enough to form a dynamically stable system.

Our detection of signals with semi-amplitudes to below 0.5 m/s demonstrates the ability of the RV method to find Earth analogs orbiting Sun-like stars. This work supports the conclusion of Dumusque et al. (2017) that the combination of Bayesian methods with the modeling of correlated noise is essential in finding small planets. In particular, the noise model framework developed by Tuomi et al. (2013b) and Feng et al. (2017b) and applied in this work plays a key role in disentangling small signals from noise.

Acknowledgements

The authors are supported by the Leverhulme Trust (RPG-2014-281) and the Science and Technology Facilities Council (ST/M001008/1). We used the ESO Science Archive Facility to collect radial velocity data sets. The authors gratefully acknowledge the HARPS-ESO team’s continuous improvement of the instrument, data and data reduction pipelines that made this work possible. Finally, the authors thank the referee, Xavier Dumusque, for his valuable comments that helped improving the manuscript significantly.

Appendix A Residual visualization