Color difference makes a difference: four planet candidates around tau Ceti

The removal of noise typically correlated in time and wavelength is one of the main challenges for using the radial velocity method to detect Earth analogues. We analyze radial velocity data of tau Ceti and find robust evidence for wavelength dependent noise. We find this noise can be modeled by a combination of moving average models and"differential radial velocities". We apply this noise model to various radial velocity data sets for tau Ceti, and find four periodic signals at 20.0, 49.3, 160 and 642 d which we interpret as planets. We identify two new signals with orbital periods of 20.0 and 49.3 d while the other two previously suspected signals around 160 and 600 d are quantified to a higher precision. The 20.0 d candidate is independently detected in KECK data. All planets detected in this work have minimum masses less than 4$M_\oplus$ with the two long period ones located around the inner and outer edges of the habitable zone, respectively. We find that the instrumental noise gives rise to a precision limit of the HARPS around 0.2 m/s. We also find correlation between the HARPS data and the central moments of the spectral line profile at around 0.5 m/s level, although these central moments may contain both noise and signals. The signals detected in this work have semi-amplitudes as low as 0.3 m/s, demonstrating the ability of the radial velocity technique to detect relatively weak signals.


Introduction
The radial-velocity (RV) technique is one of the most successful methods used to detect exoplanets. The extreme precision spectrographs developed in recent years have improved the precision of Doppler measurements down to a few meters per second. In particular, the High Accuracy Radial Velocity Planet Searcher (HARPS) spectrometer has enabled the discovery of Super-Earths due to its precision of measuring down to 1 m s −1 RV (Pepe et al. 2002;Mayor et al. 2003). However, this precision is still not high enough to detect Earth analogues in the habitable zone of nearby stars, which requires achieving 10 cm s −1 precision (Mayor et al. 2014). Moreover, efficient statistical tools and noise models are required to disentangle the signals from stellar and instrumental noise, as summarized in the RV challenge results (Dumusque et al. 2017).
The Keplerian signals in RV measurements can be diluted and distorted by stellar activity, rotation, and uneven sampling of observation times. These sources of contamination can be partly removed by various activity indicators such as the Ca II HK emission (RHK), line bisector span (BIS), and the width of the spectral lines (FWHM). However, the relation between the indicators and their RV counterparts could be very complex and is not necessarily deterministic, leading to controversial results in the validation of planetary candidates (e.g., Robertson et al. 2014;Anglada-Escudé & Tuomi 2015). This incomplete modeling of RV noise and lack of consensus on the most appropriate and efficient statistical methods are limiting the abilities of RV analysis to detect Earth analogues (see Feng et al. 2016 for details). Nevertheless, our noise modeling approach is one of the best RV modeling strategies in the field according to the analysis of the RV fitting challenge results (Dumusque et al. 2017).
Another challenge is the dependence of RV noise on wavelength or, in practice, on echelle order because RVs for HARPS are determined on an order-by-order basis. Since the jitter in RV variations depends on spectral orders (Anglada-Escudé & Butler 2012), the RV averaged over all orders would contain wavelength-dependent noise due to a lack of appropriate weighting and correcting. This motivates us to model the color dependency of the RVs. We divide the 72 spectral orders into groups, and average the RVs in each group to generate the so-called "aperture data sets" and investigate the differences between these aperture data sets-the so-called "differential RVs" (Feng et al. 2017b).
In this work, we use a combination of moving average models and the differential RVs to remove wavelength-and time-dependent noise. We apply this model to the HARPS measurements of τ Ceti, which may host a multi-planetary system according to previous analyses hereafter MT13). τ Ceti is a Sun-like star but is not as active as the Sun. There are currently more than 9000 HARPS measurements of this star, potentially enabling us to find signals with semi-amplitude as low as 0.2 m s −1 (MT13). Although MT13 have removed part of the correlated noise using moving average models, their noise modeling is probably incomplete, since the wavelength-dependent noise was not taken into account. With new data obtained by HARPS and the use of differential RVs, we reanalyze the RV variations of τ Ceti to find Keplerian signals and to attempt to verify the results of MT13. This paper is structured as follows. First, we introduce the HARPS and Keck measurements of τ Ceti, and define various data sets in Section 2. In Section 3 we describe the Markov Chain Monte Carlo (MCMC) method used to sample the posterior distribution within the Bayesian framework. Then we justify the use of differential RVs to remove wavelengthdependent noise, and select the optimal noise model for each data set in Section 4. We apply these models to find planetary candidates, and compare them with previous results in Section 5. We also investigate the cause of highly eccentric signals. In Section 6, we report the parameters of planetary candidates, and analyze the dynamical stability and habitability of these planetary candidates. Finally, we discuss and conclude in Section 7.
The main data we will use are the RVs measured by HARPS (Mayor et al. 2003) and processed by the TERRA pipeline (Anglada-Escudé & Butler 2012). The data are processed using the astrocatalog mode of TERRA whereby all barycentric corrections are recomputed using consistent ephemeris and coordinates and proper motions based on van Leeuwen (2007). This means that the calculation of barycentric earth RV does not rely on telescope header information input by the different HARPS programs from which we have used data. The TERRA algorithm also produces 72 data sets, one for each HARPS spectral order. Each of them is composed of RVs measured at a certain wavelength range. The RVs are analyzed in combination with three activity indices, including the calcium activity index (S-index or I S ), the spectral line bisector (BIS or I B ) and FWHM (or I F ) of the cross-correlation function (CCF). To increase the signal-to-noise ratio in aperture data sets, we evenly divide the 72 RV orders into groups, and average the data sets by order in each group weighted by their measurement uncertainties to form an averaged data set. For example, we divide the 72 orders into n groups, and average the data sets in each group to generate n aperture data sets, named nAPi, where i=1,K, n. The average of all orders forms the 1AP1 data set.
To remove short-term noise in the RV data sets, we define another type of RV by binning the RVs measured within one hour. We start from the beginning of an RV data set and set the beginning time as the reference time. Then we average the RVs within one hour from the reference time weighted by their uncertainties. We then define the first time point out of the one hour window as the next reference time, and average the RVs within the one hour window in the same way. Repeating this, we generate the binned version of a given RV data set. The binned version of the nAPi data set is dubbed "binnednAPi." The outliers beyond 5-σ of the RVs in 1AP1 are removed from all aperture data sets. Considering that the noise caused by stellar activity may not be properly estimated by measurement errors (including photon poisson noise and a calibration error of 30 cm s −1 ), we also weight each data set by a constant based on the sum of jitter and measurement uncertainty terms. We try different jitter levels and do not find significant changes in the periodograms of the aperture data sets. In other words, the signals in the aperture data sets are not sensitive to the weighting function used to average spectral orders. Thus we still use the measurement errors to weight data sets in the averaging process.
We define the RV differences between aperture data sets as differential RVs. We denote them by "nAPx-x′," where n is the number of divisions, and x and x′ denote different data sets in the n divisions of the 72 aperture data sets. For example, by subtracting 3AP1 from 3AP2, we obtain the 3AP2-1 data set. We will use differential RVs to remove the wavelengthdependent noise in Section 4.
Apart from the TERRA-reduced HARPS measurements of τ Ceti, we also use the HARPS data reduced by the CCF method and the RVs measured by the HIRES spectrometer on the Keck telescope (Butler et al. 2016). For HIRES we model the dependence of RV variation on the photon count and integration time. The 1AP1 and Keck data sets and their normalized activity indices are shown in Figure 1. These data together with the aperture data sets are published electronically. Around JD 2453280, the FWHM scatters greatly and the RV changes rapidly. From JD 2453280 to JD 2453285, τ Ceti was observed for asteroseismology purposes. The star was observed continuously for five days, with exposure times of 40 s (Teixeira et al. 2009). For such short exposure time and high-cadence measurements, the data is contaminated by excess noise from a periodic guiding error (Teixeira et al. 2009). So we remove the 1597 data points before JD 2453500 to form a more conservative subset named "C1AP1." The whole process of data reduction and modeling is shown in Figure 2.

Data Analysis Methods
We compare RV models in the Bayesian framework. We start from the Bayes theorem, which is where M i is a model, and  is the data, P M i  ( | ) is the posterior, and P M i  ( | ) and P M i ( ) are the evidence (or marginalized likelihood) and the prior of model M i , respectively. The denominator of Equation (1) is a normalization term. The posterior of a model is a measure of its plausibility. If no model is preferred a priori, the posterior ratio is equal to the evidence ratio, which is also called the "Bayes factor" (BF). We claim that a model is favored over another if the BF is larger than a certain value. According to the analyses of the RVs for M dwarfs by Feng et al. (2016), the BF estimated by the Bayesian information criterion (BIC), combined with a threshold of 150, avoids false positives and negatives. 6 This is equivalent to the criterion, BIC 10 D > , suggested by Kass & Raftery (1995). Although they did not test this for RV data, we have found it to work well. Thus we use the BIC to estimate the BF, and select signals using the BF threshold of 150 in this work. For a given model M and data , we need to estimate its parameters q by calculating their posterior densities according to and P M q ( | ) are the likelihood and prior distributions, respectively. The specific likelihood and prior distributions for various models will be introduced in the next section. Because the posterior distributions for RV models are always multi-dimensional and multi-modal, the prior sampling may not well resolve the narrow posterior maxima. Thus we sample the posterior using the MCMC implemented by the adaptive Metropolis-Hastings algorithm (Haario et al. 2001), which was first applied to the analysis of RV data by Tuomi & Jenkins (2012). We first launch tempered/hot chains to explore the whole parameter space and to find local posterior maxima. We then generate untempered/cold chains to explore these maxima in order to quantify the signals and to estimate the parameters "maximum a posteriori" (MAP). Our method is similar to the parallel tempering MCMC algorithm introduced by Gregory (2005).
To make the MCMC chains identify the most probable areas of the posterior, we generate hot and cold chains in the following way. First, we run cold chains to obtain the posterior density for the null hypothesis that no signals exist in the data. Second, we evenly divide the logarithm of the period range into 20 intervals, and launch four hot chains to find the posterior maximum for each. We call all of the maxima found by these chains primary signals. Third, we choose the primary signal that gives the highest likelihood, and start cold chains from the position of this signal to generate a statistically representative posterior sample. Fourth, we compare the one-planet model with the null hypothesis using the BF threshold. If the one-planet model is favored, we move on to add another primary signal into the model, and run MCMC chains to calculate the BF for model comparison. We run steps 3-4 repeatedly until an extra Keplerian signal is not favored by the data. In addition to the BF threshold, we follow Tuomi (2012) to confirm a signal if it is constrained from above and below in the marginalized posterior distribution, i.e., P P M , k  ( | ) converges to a stationary distribution.

Wavelength-dependent Noise Models
The jitter in RV variations is probably wavelength dependent due to the wavelength-dependent modulation of stellar radiation by stellar activity (Desort et al. 2007;Huélamo et al. 2008). This dependence indicates that the RVs determined through averaging RV orders are biased due to a lack of proper weighting and correction. To avoid this bias and model wavelength-dependent RVs, we introduce a new type of noise proxy called "differential RVs," defined by the RV differences between aperture data sets. Since Keplerian signals do not depend on wavelength, the differential RVs only contain wavelength-dependent noise. Actually, they are approximations of the derivatives of RV noise with respect to wavelength. They provide important information about stellar activity and instrumental noise, and thus can be used to remove the wavelength-dependent noise in Doppler measurements.
If we express the non-Keplerian part in an aperture data set with wavelength range centered at λ as t, l Y( ), the differential RV is defined by the difference between two aperture data sets, ) is the RV at time t i and wavelength j l . If the wavelength difference is small enough, D t , ) is approximately proportional to the first partial derivative of Ψ with respect to j l , Similarly, if the wavelength difference is a constant ( l D ) for all differential RVs, the second partial derivative of Ψ is In principle, we can continue the above calculations to derive higher order derivatives of Ψ from differential RVs. These derivatives can be used to approximately reconstruct the wavelength-dependent noise at a given time. For example, if the dependence of noise on time and wavelength is described by a quadratic polynomial function with time-varying parameters, i.e., the differential RVs could approximate p(t) and q(t) according to (4) and (5). Hence the discrete form of the above equation is Considering that the dependence of noise on time and wavelength is probably much more complex than the above case, we estimate the noise using the following equation: where N λ is the number of aperture divisions and thus N 1 l is the number of independent differential RVs, d m characterizes the linear dependence of RV noise on differential RV at wavelength m l , which is averaged over the differential RV, a is the acceleration caused by activity cycles or other wavelengthdependent noise 7 and b is the reference velocity. The linear dependence of RVs on activity index I k is parameterized by constant c k . All these parameters are wavelength dependent up to a certain level. This equation is equivalent to a set of independent noise models applied to different aperture data sets.
To predict the RV variation, we combine n Keplerian components with a wavelength-dependent noise component to form a basic RV model, which is , , ( ) is the RV variation caused by the kth planet, and Ŷ is the estimation of the noise component. In the Keplerian component, K k , k w , k n , and e k are the amplitude, longitude of periastron, true anomaly, and eccentricity for the kth planetary signal. According to Feng et al. (2016), the time-varying RV noise (r(t) in our case) of M dwarfs can be well modeled using a combination of white noise and the first-order moving average (MA) models. However, τ Ceti is much hotter than M dwarfs and thus the convective velocity is greater, potentially giving rise to significant granulation and stellar oscillation signals (e.g., Kjeldsen & Bedding 1995;Meunier et al. 2017). Moreover, the data is sampled with high cadence and thus probably contaminated by significant correlated noise. Thus we consider higher order moving average models to remove the RV noise, following MT13. We define the general moving average model with exponential smoothing as where w k and τ are the amplitude and timescale of the moving average while i k  -is the residual after subtracting the data by a realization of the basic model at time t i k -. Hereafter, we define "nP+MA(q)+mD" as the n-planet model with qth-order moving average and m differential RVs, which are derived from m 1 + aperture data sets (see Figure 2). The white noise model is denoted by MA(0).
For the 1AP1 data sets, the wavelength j l should be regarded as an averaged wavelength. Equation (9) is still applicable because the averaging process is a linear process, and thus the linearity of the RV model remains. Although the moving average model has an exponential term, the correlation timescale τ is not sensitive to wavelength according to our analysis. Hence we will apply Equation (9) to 1AP1 and other aperture data sets.
For a given aperture data set, the excess white noise or white jitter is taken into account in the likelihood s is the measurement noise at time t i in the jth aperture data set, which has an averaged wavelength of j l , s J j l ( ) is the jitter level, and v i j , is the observed RV at time t i in the jth aperture data set.
Following Feng et al. (2016), we adopt uniform prior distributions for most parameters except for the eccentricity and some timescale parameters. Since the planets with highly eccentric orbits are very rare, we adopt a Gaussian distribution centered at zero and with a standard deviation of 0.1 (Tuomi & Anglada-Escudé 2013). Although we set a lower limit of one day for the orbital period, we do investigate shorter period signals if the power of short period is high in the periodogram. The prior distributions of all parameters are described in Table 1.
Nearly all RV data analysis has only considered the averaged data or 1AP1 is analyzed without accounting for wavelengthdependent noise. Here, we demonstrate the wavelength dependence of RV noise by applying the model defined by Equation (10) to aperture data sets in Sections 4.2 and 4.3. In the other sections, we will apply Equation (10) to the 1AP1 and other data sets in order to identify Keplerian signals. The number of Keplerian components, MA components, and differential RVs will be chosen based on Bayesian model comparison.

Removing Wavelength-dependent Noise
To see whether the differential RVs can reduce the wavelength dependence of RV noise, we compare RV models with and without differential RVs in the Bayesian framework. We apply the RV model defined in Equations (10) and (11) to the 3AP1, 3AP2, 3AP3, and 1AP1 data sets. Specifically, we model the aperture data sets using the first-order moving average and white-noise models without and with dependence on differential RVs, which are denoted by 0P+MA(1)+0D, 0P+MA(1)+2D, 0P+MA(0)+0D, and 0P+MA(0)+2D, respectively. To visualize the differences between various noise models, we calculate the generalized periodograms for the data sets and their residuals.
The generalized periodogram is calculated by maximizing the logarithmic likelihood of a combination of sinusoidal functions and a linear trend at a sample of frequencies. Unlike the Lomb-Scargle periodogram, this general periodogram not only optimizes the amplitude and phase in the sinusoidal function but also optimizes the reference velocity and linear acceleration for each frequency. We call this periodogram the generalized Lomb-Scargle periodogram with floating trend (GLST 8 ; Feng et al. 2017a; also see Baluev 2008 andSüveges et al. 2015), a generalization of the so-called generalized Lomb-Scargle periodogram (GLS; Zechmeister et al. 2009).
Following Cumming et al. (1999), the GLST is normalized using the residual variance, and the FAP values are calculated accordingly. They FAP is not accurate for signal identification/quantification especially when the RV noise is highly correlated in time and wavelength. Thus we only show FAPs in the GLSTs of residuals where correlated noise is properly modeled and subtracted. Although there are calculations of more precise FAPs (Baluev 2008) and more sophisticated Linear Trend and Jitter Activity Indices and Differential RVs Note. The unit of c k and d m is m s −1 because the activity indices and differential RVs are normalized to zero mean and unit standard deviation before inclusion in the model. The maximum and minimum times of the RV data are denoted by t max and t min , respectively. The maximum amplitude of the RV data set with respect to the mean is denoted as v v max -|¯| . The parameter characterizing the dependence of RV on activity indices is c k , where k denotes the names of various indices. Figure 3. We aim to illustrate the existence of wavelength-dependent noise and the necessity of removing it using differential RVs. The interpretation of specific peaks in plotted periodograms is not important. Rather, our concern is the consistency between periodograms of RVs and their residuals measured at different wavelengths. Thus we show a series of different periodograms across the page. Each row of plots is given to a different wavelength range. Going down the page they are for the data sets 3AP1, 3AP2, 3AP3 (splitting the data set into three parts), and 1AP1 (regular data set averaged over all orders). The following columns show periodograms for the corresponding residuals after subtracting the model predictions for the 0P+MA(0)+0D (white noise), 0P+MA(1)+0D (moving average), 0P+MA(0)+2D (white noise with differential RVs), and 0P+MA(1)+2D (moving average with differential RVs). The logarithm BF of a model with respect to the MA(0) model are shown for the residuals of each data set after subtracting the best model prediction. The red dotted lines denote the periods at the maxima of posterior distributions. periodograms like the ones introduced by Baluev (2013Baluev ( , 2015 and by Feng et al. (2017a), we rely on Bayesian posterior samplings to identify and quantify signals.
The periodograms for the 3AP1, 3AP2, 3AP3, and 1AP1 data sets and their residuals after subtracting the model prediction are shown in Figure 3. Notably we observe great differences between periodograms for the data sets of 3AP1, 3AP2, 3AP3, and 1AP1, which indicates that wavelengthdependent noise is an important factor. On the contrary, the residuals of all RV data sets have similar periodograms after subtracting the best predictions of 0P+MA(0)+2D or 0P+MA (1)+2D from the data, indicating the essential role of differential RVs in removing wavelength-dependent noise. In the right two columns of panels in Figure 3, we see consistent periodograms for residuals despite being calculated for different noise models and aperture data sets.
For the 1AP1 data set, the differential RVs do not improve the BF as much as the moving average model. On the other hand, the increase of logarithmic BF by including both MA(1) and differential RVs in the model is approximately equal to the sum of those by including them separately (see the left three columns). This indicates that the wavelength-dependent noise and time-correlated noise are independent and thus should be modeled with independent noise components. The inclusion of differential RVs improves the fitting more for the 3AP1 and 3AP3 data sets than it does for the 3AP2 and 1AP1 data sets. This means that the RVs measured at the middle of the wavelength range contain less wavelength-dependent noise than those measured at the blue and red ends of the range do. The 1AP1 data set contains less correlated noise probably due to a partial removal of wavelength-dependent noise by the averaging of all spectral orders. But this noise is still significant enough to contaminate the periodogram and thus result in detections of noise-induced signals.
In summary, for both the 1AP1 and the other data sets, the differential RVs are able to remove wavelength-dependent noise and help avoiding false positives. This role of differential RVs is rather different from that of red noise models, which are good at removing time-correlated but not wavelength-correlated noise. We also conduct similar analysis for the binned data sets, and find similar results, although the wavelength-dependent noise is greatly reduced by the binning process. Considering that τ Ceti is a quiet main sequence star, our analysis is probably representative. The wavelength dependence of RV noise is probably stronger for more active stars such as Alpha Centauri A/B (Dumusque et al. 2012). For them, differential RVs would play a key role in detecting exoplanets consistently.
In the following sections, we will focus our analysis on the 1AP1 data set because it has higher signal to noise than aperture data sets, and has residuals similar to those of aperture data sets. We further study the effect of binning by comparing the periodograms for the binned1AP1 and 1AP1 data sets and for their residuals in Figure 4. We find that periodograms are very different both for the RV data sets and for their residuals. The binning of data over a one hour time span has removed features caused by noise and signals altogether. This makes the binned version unreliable for detecting signals.

Modeling Instrumental Noise
Like stellar activity, instrumental noise may also be colorful because of the instrument's potential nonlinear response to wavelength. To model the instrumental noise in the data, we generate the calibration data sets by combining the RV aperture data sets measured by HARPS and reduced by the TERRA algorithm for 172 stars that have been reduced with the TERRA reduction and represent the most frequently observed targets (excluding τ Ceti). Most of these RVs were measured within the HARPS-Upgrade GTO program (Mayor et al. 2011;Pepe et al. 2011). We derive aperture data sets and differential RVs from these HARPS measurements, and combine them. Specifically, we remove the (differential) RVs that have absolute values larger than 20 m s −1 or deviate from the mean more than 5s before combining them. For each epoch in each aperture data set for τ Ceti, we average the calibration data points that have nearest epochs by weighting them according to their measurement errors. We further remove the outliers that deviate from the mean more than 3s. There are also epochs where no RVs of other stars are available; we assign the RVs measured at nearby epochs. We use these calibration data sets as proxies (like activity indices) to remove instrumental noise. For example, we can use a linear combination of the calibration data sets for C3AP2-1 and C3AP3-2 to model the instrumental noise in the C1AP1 data set. Hereafter, we use cCnAPi (cnAPi) and cCnAPi-j (cnAPi-j) to denote the calibration data sets for CnAPi (nAPi) and CnAPi-j (nAPi-j).
Through comparing various combinations of calibration data sets for all RV data sets, we find that the calibration data sets are appropriate proxies to reduce the instrumental noise in the data. Although the cC1AP1 data set is influenced by the instrument in the same way as the C1AP1 data set, the former is "contaminated" by Keplerian signals from other stellar systems though the contamination is reduced by removing outliers and averaging the data. On the contrary, the differential calibration data sets do not contain Keplerian signals, and thus are more appropriate for removing instrumental noise. We find that the dependence of RV on calibration data sets is not consistent with zero. The HARPS measurements are biased at least by 0.244 0.035 0.083 -+ m s −1 , determined by the square root of the sum of MAP estimations of linear coefficients for all cC9AP differential data sets. This bias is not reduced by including Keplerian components into the model, suggesting its instrumental origin. The real bias caused by the instrument could be higher because the differential calibration data sets only account for the wavelength-dependent instrumental noise. Thus the HARPS data would not be reliable for detecting signals below 0.2 m s −1 if the instrumental noise is not properly removed.
However, not all of the cC9AP differential data sets are useful, and the application of all of them only slightly improves the BF with respect to the cC3AP differential data sets. For the C3AP sets, only cC3AP2-1 is strongly correlated to the RV data, indicating higher instrumental noise in blue spectral orders than in red orders. In Figure 5, we show the cC3AP2-1 data set for C1AP1, and the corresponding posterior distribution for the model. We see a correlation between cC3AP2-1 and C1AP1 up to 0.15 m s −1 . Hence we will use the cC3AP2-1 data set to remove wavelength-dependent instrumental noise in the HARPS data. Hereafter we include cC3AP2-1 together with the S-index, BIS, and FWHM indices linearly in the model (see Equations (10) and (8)) to analyze the C1AP1 data set. Similarly, we use c3AP2-1 to remove the instrumental noise in the 1AP1 and CCF data sets. We select the optimal number of MA components and differential RVs in the following subsection.

Choosing the Optimal Noise Model
Within the Bayesian framework, we compare noise models with various differential RVs and MA components to determine the optimal noise model for the C1AP1 data set. We report the BFs with respect to the MA(0) model (without differential RVs or any Keplerian component) in Table 2. The optimal number of differential RVs and MA components are selected based on the BIC-estimated BF threshold of 150 (Feng et al. 2016). In other words, the optimal model is the most complex model which gives a BF at least 150 times higher than all simpler models. This is equivalent to an increase of 5 of logarithmic BF in Table 2.
According to this criterion, the 0P+MA(4)+8D model is the noise model favored by the 1AP1 and C1AP1 data sets. Including more differential RVs may not reduce the noise; rather, the noise in the differential RVs may reduce the significance of true signals. To explain this we compare the periodograms for the residuals of the 1AP1 data set after subtracting the 0P+MA(4)+8D and 0P+MA(4)+71D models in Figure 6. Guided by the power difference between FAP thresholds in the periodograms, we see that the noise in the residuals increases if we include more differential RVs. Moreover, complex noise models may interpret signals as noise due to their flexibility (Feng et al. 2016 and MT13). Based on the above considerations, we adopt the 0P+MA(4)+8D model to remove RV noise. We also find at most weak correlation between noise proxies (see Appendix A), supporting the necessity of using all of them in noise modeling. This Goldilocks model is different from the one devised by Feng et al. (2016) for M dwarfs because the target and data sets in this work are different. But results both in this work and in Feng et al. (2016) suggest the importance of finding an appropriate noise model for each specific RV data set.
Following this approach, we further find the Goldilocks noise models for the CCF data set to be 0P+MA(5)+8D. The optimal noise model for the Keck data set is 0P+MA(1)+0D. In the 0P+MA(4)+8D model; there are two free parameters for trend, one parameter for jitter, four parameters for activity indices and c3AP2-1, five parameters for the MA model, and eight parameters for differential RVs. Thus there are 20 and 21 free parameters in the 0P+MA(4)+8D and 0P+MA(5)+8D models, respectively. Since the differential RVs can weight aperture data sets a posteriori and thus remove wavelengthdependent noise (see Appendix B), it is unnecessary to analyze aperture data sets separately. Hence we only analyze the averaged data sets, C1AP1, 1AP1, and CCF, to identify signals.

Performance of the Goldilocks Noise Model
Based on visual investigation of the 1AP1 data set, we see a rapid increase in RV around epoch JD 2453282 and decrease around epoch JD 2456190, as shown by black points in  . We see a steady increase of RV by 5 m s −1 and a rapid decrease by 10 m s −1 around the above epochs. These variations are huge compared with the sub-meter semiamplitudes of the signals, which we will report in the following section. Any model that does not fit these two features would have low likelihood, and thus not be favored by the data.
We use the Goldilocks noise model, 0P+MA(4)+8D, to model the noise in the 1AP1 data set and show the binned residual after subtracting the best-fitting model as red points in Figure 7. We see that the noise model significantly reduce the noise-induced RV variation. However, we still see weak intranight noise in the epochs around JD 2453283, probably caused by the guiding error, as mentioned in Section 2. Hence the C1AP1 data set is probably a more conservative choice for signal detection. According to our analysis, the Keck data does not help to constrain the signals detected in HARPS data sets. To be conservative, we will analyze C1AP1 and CCF to identify signals, and use Keck and 1AP1 for the sensitivity and consistency test in the following sections.

Keplerian Signals
In the above section, we obtain the optimal noise model for the TERRA-reduced HARPS data set. Since the spectral orders of the CCF data is not available, we apply the TERRA differential RVs to model the RV noise in the CCF data set. 9 From now on, we only use the C1AP1 data set to identify signals and use 1AP1 and CCF to test the consistency of signals.

Primary Signals
As mentioned in Section 3, we run hot chains to find primary signals for later investigations. We divide the period range into 20 chunks and run a hot chain to find the local posterior maxima for each. We convert the tempered posterior to the posterior for each chain, and combine all posteriors to approximate the posterior distribution for the 1P+MA(4)+8D model. We also start cold chains from the local posterior maxima of hot chains, and combine the posteriors drawn by cold chains. To compare with signals detected in the Keplerian solution, we also fit sinusoidal functions to the data to obtain circular solutions. In Figure 8, we show the distributions of logarithmic BF (estimated by BIC at a given period) over period for the Keplerian and circular solutions for the 1AP1 data set. We see that the strongest signals are around 160, 600, and 1000 days. The 600 and 1000 day signals are probably annual aliases of each other. The 114 and 318 day signals are annual aliases of 160 days, while 226 days is an annual alias of 600 days. The signals at periods of about 20, 49, 160, 600, and 1000 days and some of their annual aliases are also significant in the Keplerian solution.
To demonstrate the uniqueness of the primary signals, we show the posterior distribution for two period ranges where the signals around 20 and 49 days are identified in Figure 9. We see that the two signals are unique and well identified by the hot chains in the corresponding period intervals. Therefore it is reliable to use the parameters of these signals as initial conditions for cold chains to constrain signals simultaneously.
In Figure 10, we show the logarithmic BF distributions of the samples drawn by cold chains for the C1AP1 and CCF data sets. For both solutions and data sets, the signals around 160, 600, and 1000 days are most significant. The 600 day signal is more significant than the 1000 day signal for both Keplerian and circular solutions for both data sets. According to our analysis, the signals at periods of 20, 49, 160, 600, and 1000 days can also be found in the 1AP1 data set and combinations of HARPS and Keck data sets. However, the comparison of the significance of primary signals is biased because the one-planet model may interprets the variations caused by multiple signals as one Keplerian signal. Thus we will constrain these primary signals simultaneously, and compare the results for various data sets in the following section.

Comparing Signals Detected in Different Data Sets
Using the numerical method described in Section 3, we obtain posterior samples from the MCMC sampling for RV models with various numbers of Keplerian components. We Figure 5. The normalized cC3AP2-1 calibration data set (left), and the posterior densities of the linear dependence of the C1AP1 data set on it (parameterized by c c ) for the 0P+MA(4)+8D model (right). The calibration data sets are normalized to zero mean and unit standard deviation before inclusion in the model. Thus c c is in m s −1 , characterizing the level of instrumental bias. The mode, mean, standard deviation, and higher order moments of the posterior density are also shown. 9 Although these differential RVs are produced by TERRA, they should be able to remove wavelength-dependent noise in the CCF data since the two data sets are rather similar, leading to at most a second-order difference in the differential RVs produced by TERRA and CCF. Even though they produce rather different differential RVs, we can still regard TERRA differential RVs as a type of activity indices and consider their linear correlation with the CCF RVs. select signals according to the signal detection criteria within the Bayesian framework. The parameters of signals identified in CCF and C1AP1 are shown in Table 3, and the BFs for the Keplerian and circular solutions for the C1AP1 set are shown in Table 4. The signals detected in the Keck data set are not shown because only a signal at 20 days is identified. The strongest signals in Keplerian solutions for both data sets have orbital periods around 1000 days and 600 days. Since they are annual aliases, we report the results for both in the table. We find that the eccentricities of these two signals are high for CCF but are below 0.2 for C1AP1. This indicates that the high eccentricity is probably caused by the noise in the observation epochs before JD 2453500 (see Figure 1), which are included in CCF but not in C1AP1.
Highly eccentric and low-mass exoplanets are very rare (Tuomi & Anglada-Escudé 2013), which is evident from the distribution over the mass and eccentricity of all exoplanets detected through the RV technique in Figure 11. Although there are many biases in a plot like this and particularly at low masses and high eccentricities, we do see a lack of planets with low mass, long period, and high eccentricity like the 1000 day signal we have identified in this work. Thus the higheccentricity solutions reported in Table 3 are probably not caused completely by planets orbiting the star but are superpositions of planets and activity. The activity and signals are not completely disentangled probably because of incomplete noise modeling or inefficiency of the MCMC sampling, which makes the chain unable to jump out of local maxima to find extra low eccentric signals. We investigate the former by adding more MA components to the noise model, but fail to reduce the eccentricity. Then we investigate the latter by finding signals with zero eccentricity, which makes the MCMC chains achieve convergence more efficiently due to reduction in dimensionality. Based on the Bayesian model comparison, the circular solutions identify more signals than the Keplerian solutions because the sinusoidal function has fewer parameters than the Keplerian function and thus is less penalized by the BIC. We find consistent circular solutions of five signals for all data sets. The 92 and 102 day signals are probably an alias pair since the subtraction of one from the data would weaken the other. The 92 day signal may be genuine because it is identified in circular solutions as well as the Keplerian solution for CCF. Notably the posterior of the 92 day signal is rather low in the Keplerian and circular solutions because it does not pass the BF threshold as a primary signal, as shown in Figure 10.
For both of the circular and Keplerian solutions, we find that the 600 day signal is favored over the 1000 day signal for all data sets. This also leads to a slightly higher amplitude of the 600 day signal with respect to the 1000 day signal (see Table 3). Moreover, the 160 day signal is very eccentric (e=0.96) in the 1000 day solution for C1AP1. Thus the longperiod signal probably has a period around 600 days, if it is caused by a planet. This signal could be confirmed by further observations of the inner edge of the debris disk of τ Ceti by the Atacama Large Millimeter/sub-millimeter Array (ALMA; MacGregor et al. 2016). If the inner edge is found to be beyond 1.5 au, this planet must exist to clear the inner disk. As for the other signals, we regard the 20 day signal as a genuine planetary candidate since it is also identified independently in the Keck data set. The signal around 160 days is Keplerian because it is the most significant signal in the circular solutions and is consistently identified in all HARPS data sets. The signal at a period of 49 days is Keplerian because it is identified in all data sets. There could be a signal at a period of about 92 days or 102 days. But the confirmation of this requires further observations and analysis. These results are also valid for the 1AP1 data set and combinations of HARPS and Keck data sets according to our analysis.

Curse of Eccentricity
Although the eccentricity of the 600 day signal is low for the C1AP1 data set, the other signals have eccentricities that are not consistent with zero. We also add extra Keplerian signals to the Keplerian solutions for the C1AP1 and CCF data sets, but the high eccentricity of signals persists. To decide whether the signals are time dependent, we evenly divide the time span of the 1AP1 data set into three chunks, and obtain Keplerian solutions for each. We find that the jitter level for the first chunk is about 0.3 m s −1 higher than the other two chunks. The jitter may significantly dilute the signal in the first chunk because they have comparable semi-amplitude, which is around 1 m s −1 . The 1000 day signal is identified in the second and third chunks, although it is rather eccentric in the third chunk. The 600 day signal is identified in the third chunk. The 49 and 160 day signals are identified in the third chunk. The 20 day signal is found in all chunks while the 14 day signal is only significant in the first chunk. This is probably the reason why MT13 has identified the 14 day signal rather than the 20 day one based on analyses of early data. This apparent inconsistency of signals is unsurprising because they are signals with semi-amplitudes of at most 0.5 m s −1 and require large samples to confirm. Although chosen based on analyses of the whole HARPS data, the 1P+MA(4)+8D model is applied here to different chunks, which probably lead to false negatives.
We investigate this by dividing the 1AP1 data into three chunks, applying the 0P+MA(2)+2D noise model to each, and constraining signals using the whole data set. In the periodogram for the residual after subtracting the noise model prediction, we are able to identify all signals despite the noise Note. The relative maximum logarithm likelihood for each model with respect to the white noise model (i.e., MA (0)) is shown. The BF of the Goldilocks noise model is shown in boldface.
parameters varying with chunks. This suggests that the signals we have identified are probably not caused by activity cycles.

Instrumental Bias of HARPS
The above investigations of the cause of high eccentricity are concerned with the removal of activity-induced noise from the data. In this section, we will study the excess instrumental noise caused by the instability of HARPS in the sub-m/s regime. We model this noise using the moments of the HARPS line profiles, which probably reflect the flux loss caused by instrumental effects such as guiding errors (Berdiñas et al. 2016). Following Berdiñas et al. (2016), the moments are calculated in the following steps: 1. The fluxes of all spectra are scaled such that the relative flux in each echelle order is the same. 2. The spectrum recorded in all echelle orders is deconvolved to obtain a mean line profile F i (t) in absorption. After subtracting the residual continuum w(t), the line profile is inverted and normalized to be a flux distribution represented by N flux values f i (t) and RVs v i (t). 3. The moments are calculated according to the following equations, Although the central moments attempt to minimize the mean velocity contribution by subtracting M 1 , M 1 is not perfect at removing Keplerian velocities. So we only use central moments to investigate the instrumental noise rather than to quantify signals before performing any robust tests on the connection between moments and RVs. We linearly combine M 0 , M 3c , M 4c , and M 5c , 10 activity indicators, calibration data sets, and differential RVs in the model, and find strong correlations between RVs and these moments. For example, we use this new noise model to constrain the 160 day signal, which is identified in the C1AP1 data set. We find that the linear dependence of RVs on the normalized M 0 , M 3c , M4c, and M 5c are around −0.17±0.02 m s −1 , −0.62±0.02 m s −1 , 0.16±0.02 m s −1 , and −0.28±0.02 m s −1 , respectively. 11 In particular, these moments increase the logarithm BF by more than 1000 for all data sets. We show this strong correlation between RVs and moments for the C1AP1 data set in Figure 12. We observe strong anti-correlations between RVs and odd-order moments that have terms proportional to M 1 -. This indicates that the subtraction of M 1 from velocities (see Equation (14)) may not completely remove the Keplerian components, leading to a considerable correlation between RVs and central moments. We apply the four moments to all data sets by including them into the model linearly. We find that the strong signals shown in Table 3 are still significant with this noise model.
To explore the correlation between the moments and calibration data sets, we show their time series around JD 2453283 and JD 2456195 as done for Figure 7. We see strong correlation between RV and various proxies over timescales The rapid increase and decrease during these two epochs are clearly seen in other proxies. Thus the instrumental bias is probably the main reason why we find high eccentricity in Keplerian solutions even though the application of the noise model is able to significantly reduce such noise (see Figure 7). However, there are also strong variations in moments within one observation night, which do not appear in the RVs. This short-term variation may be caused by the change of τ Ceti's altitude, which modulates the intensity of the atmospheric absorption of the stellar light. Further studies are required to understand these variations before using the moments as noise proxies. In Figure 13, we observe that the c3AP2-1 data set decreases around JD 2456195 as the RVs do. Since no RVs of the 172 stars were measured or published around epoch JD 2453283, we assign the calibration data points at nearby epochs to them.
Although the moments correlate with RVs up to 0.5 m s −1 , the inclusion of them in the model may have both positive and negative effects on the identification and quantification of signals because of the possibly incomplete subtraction of Keplerian variation from moments. Considering these problems, we have only used them to test the sensitivity of signalto-noise models.

Four-planet System
According to the analyses in Section 5, we conclude that the signals with periods of 20, 49, 160, and 600 days are genuine Keplerian candidates. To further confirm them as Keplerian candidates, we check whether these signals could be caused by activity. We show the periodograms of the C1AP1 and Keck data sets, activity indices, differential RVs, and the signals in Figure 14. Although there are some peaks around 600 and 1000 days in the periodograms of observation times of the Keck data, 9AP6-5 and 9AP2-1, they are either minor peaks or only strong in one or two noise proxies. Notably the 1000 day signal is rather strong in the periodogram of 9AP6-5 despite a considerable probability of random overlap between the signals in the data and the powers in noise proxies. The 160 day signal is strong in the S-index but does not appear in other proxies. The periodograms both for C1AP1 and for Keck show strong peaks around 20 days, strongly suggesting it as a genuine Keplerian signal. Other signals are not found to be strong in the Keck periodogram probably due to the large uncertainty and small sample of Keck measurements. These signals have mean semi-amplitudes below 0.5 m s −1 , 12 which is beyond the detection ability of Keck. Since the 600 day signal has low eccentricity in the Keplerian solution for the C1AP1 data set, we show the parameters of the five signals identified in it in Table 5. The boundary of uncertainty intervals, 1% q and 99% q , are determined at the cumulative marginalized posterior probability of 1% and 99%. Specifically, they are determined by , where min q is the minimum value of parameter θ. 13 The uncertainties of parameters are probably larger than the reported values because the mean values of parameters slightly depend on the choice of data set as shown in Table 3. In Table 5, we observe that all five planets have minimum masses less than 5 M Å . In particular, τ Ceti g and h have minimum masses comparable with the Earth. τ Ceti h, e, and f have considerable eccentricities. The causes could be the instrumental noise we have investigated in Section 5 and the bias in the estimation of eccentricities since eccentricity is positive definite (Zakamska et al. 2011). The planetary masses would be double the minimum masses if the best-fitting inclination for the debris disk of τ Ceti (around 30°) is equal to the inclination of the planetary system (Lawler et al. 2014). Our detection of these small planets at such long orbits demonstrates the potential power of the RV technique in finding Earth analogues. The RV variations caused by these planets have mean semi-amplitude as low as 0.3 m s −1 , which is close to the 0.1 m s −1 limit required for detecting Earth analogues around solar analogues.
If we regard the measurements in a 15 minute time bin as independent observations (e.g., Mayor et al. 2003 andO'Toole et al. 2008), there would be 662 observations for the C1AP1 data set. The K N K RV N rms obs =´ratio for signals with K=0.3 m s −1 is 7.3, where RV rms is the standard deviation of the RVs after removing the best fitted trend and the correlation with noise proxies. This is close to the detection threshold of 7.5 based on the RV challenge results (Dumusque et al. 2017). Considering that our team has detected signals with K N 5 = without announcing false positives in the RV challenge, our detection of small signals is reliable and is consistent with previous analyses.
We also show the phase-folded data and model prediction for all signals in Figure 15. We see that the data strongly support the existence of four signals.To see whether there are extra signals in the data, we subtract the best-fitting model (the 600 days Keplerian solution shown in Table 3) from the C1AP1 data, and show the GLST of the residual in Figure 16. We do not see any strong signals, supporting the fact that the six-planet model is not favored by the data.
Apart from the HARPS data set, we also combine the RVs measured by the Anglo Australian Telescope (AAT; the data is available in MT13) and Keck with HARPS. We find results consistent with the five-planet solution. Thus the weak wobbles found in HARPS data are consistent with other data despite not being identified independently due to the relatively lower RV precision of other instruments.

Comparison with Previous Studies
The planetary candidates we have identified partially overlap with the ones found by MT13. We find two new planetary candidates with periods around 20 and 49 days, but fail to confirm the signals around 14 and 35 days. Although there is evidence for the existence of the signal at around 92 days, we cannot confirm it as a Keplerian candidate because it cannot be consistently identified in all data sets and solutions. The signal around 14 days becomes weak when we subtract the 20 day signal from the data. But the opposite is not true, suggesting a non-Keplerian origin of the 14 day signal. Nevertheless, the 14 day signal does exist in some Keplerian and circular solutions after accounting for the 20 day signal. In addition, by evenly dividing the data into three chunks, we find that the 14 days is significant in the first chunk while the 20 days is significant in the second and third chunks. Since MT13 only analyzed the first two chunks, it is expected that they find the 14 day signal to be favored. However, as we have mentioned in Section 5.3, the first chunk is problematic because of its high jitter level and the abnormal variation of FWHM before JD 24503280. Moreover, this chunk contains high-cadence asteroseismology data which is influenced by guiding errors. Considering these factors, the 14 day signal is probably an activity-induced signal rather than an alias of 20 days.
The 35 day signal reported by MT13 or the 32 day signal identified in this work is very weak in the C1AP1 data set (see Figure 8). Moreover, the 35 day signal identified by MT13 is close to the rotation period of 34 days determined by Baliunas et al. (1996), although this rotation period is not significant in the periodograms for activity indices and for differential RVs (see Figure 14). We also find that the alias of the 600 day signal, located at around 1000 days, fits the data as well as the 600 day signal if high eccentricity is allowed (see Table 4).
The semi-amplitudes of signals identified in this work are from 0.4 to 0.5 m s −1 . These values are lower than the values reported by MT13, which are from 0.58 to 0.75 m s −1 . In the analysis of MT13, the activity-induced variation is probably misinterpreted as signals because of insufficient noise modeling. Moreover, the periods of signals in this work are different from the values reported by MT13. This difference together  with the other differences are probably caused by the following factors.
First, the HARPS RV sample used in this work is double the size of the sample (4398 RV points) used by MT13. Although MT13 have also used the Keck and AAT samples, both of them have larger uncertainties and have less than 1000 measurements. The updated HARPS data set allows us to find weaker signals, and quantify signals more accurately. In addition, the early HARPS data are found to be more noisy than recent ones according to our analyses. This could also be the reason why MT13 have identified somewhat different signals.
Second, we use the differential RVs to remove wavelengthdependent noise in the data. Although the aperture data sets are not available for the Keck data, the signal is mainly constrained by the HARPS data (see Figure 15) where the noise is properly modeled by differential RVs (see Figure 3). In addition, we test different noise models and select the best one for each data set to avoid false negatives and positives.
Third, we estimate the BF using the BIC, which penalizes complex models more strongly than the Akaike information criterion (AIC) and the truncated posterior mixture (TPM) used in MT13. According to the comparison of various BF estimators by Feng et al. (2016), the BIC is more conservative in estimating the BF than the AIC and TPM estimators because it assumes that all parameters are equally free. However, this might not be appropriate due to the nonlinear combination of Keplerian parameters, leading to the unequal ability of parameters to improve the fitting. For example, the likelihood is not as sensitive to the eccentricity as it is to the period. This can be seen in Table 5 where the relative uncertainty of eccentricity is much larger than that of the period. As there is no universal method to calculate the BF, and we regard the BIC as a practical and conservative diagnostic tool to confirm signals.

Dynamical Stability and Habitability Analysis
We now assume that the detected signals correspond to planets orbiting τ Ceti. We analyze the dynamical stability of the planetary candidates using the threshold of Lagrange stability introduced by Barnes & Greenberg (2006). We calculate the instability boundary in the phase space of eccentricity and semimajor axis using the MAP estimation of the parameters of candidate planets, and show the results in Figure 17. Although the orbits of candidates are eccentric, the system is dynamically stable. As we have shown in Section 5.4, the high eccentricity is probably caused by instrumental noise. The eccentricity of the planetary candidates might be overestimated because eccentricity is definitely positive (Zakamska et al. 2011). Considering the overestimation of eccentricity, the stable region (or the gap between shaded regions) could be larger than that in the figure. On the other hand, the stable region would shrink if the true masses of these candidates are much larger than their minimum masses. But the area of the shaded regions are more sensitive to eccentricity than to the mass according to Equation (2) in Barnes & Greenberg (2006). Therefore we consider the Lagrange stability shown in Figure 17 to be a conservative illustration of the dynamical stability of the planetary system.
The four-planet system is dynamically packed, which is probably due to the effect of planet-planet scattering (Raymond et al. 2009). Our result is consistent with an inner edge as low as 1.6 au determined by ALMA observations  show that the previously reported five-planet system is stable over long timescales (Lawler et al. 2014). This indicates that the four-planet system reported in this work may also be stable for a long time since both systems are tightly packed and the outer planets have similar semimajor axes. We analyze the habitability of the planets by calculating the habitable zone according to the method introduced by Kopparapu et al. (2014). We adopt the effective temperature and luminosity of the τ Ceti measured by Santos et al. (2004) and Pijpers (2003), which are 5344 K and 0.52 L  , respectively. We find, for a one-Earth-mass planet, the conservative habitable zone of τ Ceti is from 0.70 to 1.26 au while the habitable zone for a five-Earth-mass planet is from 0.68 to 1.26 au. Thus none of the reported planets reside in this conservative zone. If one instead takes the recent Venus and early Mars limits of the habitable zone as 0.55 and 1.32 au then τ Ceti e and f are close to the inner and upper boundary of the optimistic habitable zone. However, τ Ceti e and f are probably not dynamically habitable. The ALMA observations of τ Ceti suggests the existence of a broad debris disk (MacGregor et al. 2016). The large amount of asteroids and comets in the disk may be strongly perturbed by planets, leading to an impact rate 10 times higher than that on the Earth (Greaves et al. 2004). On the other hand, the impact rate on these planets may depend not only on the mass of the scattered disc but also on the structure of the planetary system and the interstellar environment (Feng & Bailer-Jones 2014). Thus, comprehensive studies on the dynamics of the τ Ceti system are important for understanding its habitability.

Discussions and Conclusions
We analyze the RV data of τ Ceti in a Bayesian framework. We find strong dependence of the RV noise on wavelengths or orders. This colorful noise cannot be removed by averaging all spectral orders, and would cause inconsistent results if not removed. To model this noise, we introduce a wavelengthdependent noise model by linearly combining moving average models with differential RVs. We apply this model to aperture data sets, and find that the differential RVs efficiently remove wavelength-dependent noise, reducing the false positive and negative rate. We also find that binning RVs over time would change both the signals and noise in the data due to the lack of an appropriate weighting function. Therefore, we propose a combination of moving average models and differential RVs to model the time and wavelength-dependent noise in unbinned RV data sets.
Our findings challenge the traditional methods that only account for wavelength-independent noise. This means that the planetary candidates with semi-amplitudes around 1 m s −1 found in previous studies may require further confirmation by applying differential RVs to noise modeling. New data reduction algorithms are needed to extract aperture data sets from Doppler measurements for various instruments. Another application of differential RVs is concerned with the diagnostics of stellar activity. Aperture data sets provide abundant information for the study of stellar physics. The spectral powers of differential RVs are wavelength dependent (as shown in Figure 15), suggesting the colorfulness of stellar activity. Further investigations are necessary to confirm differential RVs as signatures of the spectral and photometric variations of stellar surfaces. On the other hand, our work focuses only on the analysis of the large RV data for τ Ceti. Further studies on the connection between differential RVs and stellar activity for various types of stars and instruments are needed to extend the wavelength-dependent noise model proposed in this work.
We apply the noise model to various RV data sets and get Keplerian and circular solutions. Both solutions consistently identify the signals around 600, 160, 20, and 49 days for all data sets. For the Keplerian solutions, the strongest signal has a period of either around 600 days or around 1000 days, which are aliases of each other. The 600 day signal is probably a genuine signal Note. The signals are identified in the Keplerian and circular solutions for the C1AP1 data set. The period of the long-period signal is shown in brackets. Note that the logarithm BF threshold is 5. Figure 11. Distribution of orbital periods and masses of the candidate planets detected in the C1AP1 data set (triangles) and all exoplanets detected by the radial-velocity technique (small dots). To compare the τ Ceti candidate planets with solar system planets, the parameters of the Venus, Earth, Jupiter, and Saturn are shown in open circles from left to right. The eccentricity is encoded in the color of the points. The data are collected from the Extrasolar Planets Encyclopaedia (http://exoplanet.eu).
because it is favored over 1000 days for the circular solutions. Moreover, the eccentricity of the former is low for some data sets, while the latter is highly eccentric for all data sets. Nevertheless, significantly eccentric signals exist in all Keplerian solutions. The relevant results and data are available at http://star-www.herts. ac.uk/~ffeng/HD10700_supplementary/tau_ceti_results/. Figure 14. GLSTs of the data, activity indices of Keck (top panels), C1AP1 (middle and bottom panels), and differential RVs (bottom panels). The 600 and 1000 day signals and the signals constrained in combination with the 600 day signal for the C1AP1 data set (see the fifth column of Table 3) are shown by the red dotted lines. We have truncated the periodograms at a period of 10 days to optimize visualization. By dividing the data into chunks, we see the sensitivity of data to signals varying with time. But all signals can be found if we fit the same Keplerian component and independent noise components to data chunks. Moreover, we have also studied the instrumental noise in HARPS, and find that the HARPS data is biased at the level of at least 0.2 m s −1 . The central moments of the spectral line profile are also found to be strongly correlated to the data at the level of at least 0.5 m s −1 . Future studies are required to confirm whether or not the correlation is caused by the Keplerian variation left in the central moments due to the dependence of RVs to spectral flux. The signals we have identified are robust to the inclusion of central moments in the noise model.  . GLSTs for the C1AP1 data set after subtracting the best-fitting 0P+MA(4)+8D (left) and 5P+MA(4)+8D (right) model. The signals are denoted by the red dotted lines. Since the correlated noise is subtracted, the FAPs of 0.1, 1, and 0.001 are shown by dashed lines as a metric to measure the significance of signals.
The planetary candidates found in this work partially overlap with those detected by MT13. With more HARPS measurements and wavelength-dependent noise models, we improve the estimation of the 600 and 160 signals identified by MT13. We also find two new planets together with the 1000 days alias of the 600 day signal but fail to confirm the other two planetary candidates proposed by MT13. Since the signal at a period of 92 days cannot be consistently identified in all data sets, we do not confirm it as a planetary candidate. These differences suggest that we should be very cautious about the interpretation of sub-m/s signals detected in large and high-cadence data sets. For example, the signals identified around α Centauri B (Dumusque et al. 2012) are found to arise from the window function (Rajpaul et al. 2015).
Our analysis of the dynamical stability of the four-planet system suggests that τ Ceti h, e, and f should have eccentricities lower than their values determined in this work if they are genuine planets. Although τ Ceti e and f are located close to the boundaries of the optimistic habitable zone, their habitability might be strongly reduced by the bombardments of objects from the massive scattered disc. The determination of the inner boundary of this scattered disc may help to confirm the existence of the 600 days candidate. If the scattered disc and the planetary system are coplanar, the true masses of these planetary candidates are double the minimum masses reported here (Lawler et al. 2014;MacGregor et al. 2016).
Our detection of Keplerian signals as low as 0.4 m s −1 demonstrates the unique role of the RV technique in detecting Earth analogues around Sun-like stars. Unlike the transit technique, the RV method does not require the occurrence of transits, and can be useful to detect planets around all bright stars. With the development of high-precision spectrometers, advanced statistical methods and noise models, it is feasible to find Earth analogues in the coming decade. Since the weak wobbles caused by Earth analogues are comparable to the measurement uncertainties and jitter levels, it is crucial to build appropriate noise models to avoid both false positives and negatives. Our work on modeling wavelength-dependent noise is intended as a step toward a noise model framework where physics-motivated models and stochastic processes are properly combined to remove RV noise.
F.F., M.T., and H.J. 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 teams' continuous improvement of the instrument, data, and data processing that made this work possible. Finally, the authors would like to thank the anonymous referees for their valuable comments that enabled considerable improvements of the manuscript.

Appendix A Correlation between Noise Proxies
We test the correlation between the linear coefficients of noise proxies during the fitting. We randomly select 1000 samples from the MCMC chain for the 4P+MA(4)+8D Keplerian solution (with the 600 day signal identified) for the C1AP1 data set and show the correlations between the linear coefficients of various noise proxies in Figure 18. Most noise proxies are independent of one another. There are weak correlations between differential RVs derived from nearby spectral orders because of much common noise shared by them. Nevertheless, such a correlation may not cause the so-called "collinearity" problem, considering that the MCMC chains show no evidence for non-convergence. Moreover, the application of noise proxies in fitting is seen in the analysis of Kepler light curves (e.g., Foreman-Mackey et al. 2015).

Appendix B Wavelength dependence of Model Parameters
To explore the wavelength dependence of noise model parameters, we apply the 0P+MA(4)+8D model to 9AP aperture data sets. We show the dependence of the MAP estimations of parameters and uncertainties on the mean wavelength of aperture data sets in Figure 19. We see strong wavelength dependence of jitter s J and reference velocity b. Thus a simple average of all spectral orders would introduce wavelength-dependent noise without a proper weighting function. We also see weak wavelength dependence of the trend a and FWHM c f . Notably, we observe strong wavelength dependence of differential RVs. For a given aperture data set, the coefficients for the differential RVs with shorter and longer wavelengths always have opposite signs because they are adjusted to fit the aperture data set. For example, the coefficients for 9AP2-1, 9AP3-2, 9AP4-3, 9AP5-4, and 9AP6-5 are positive while the other differential RVs are negative in the analysis of 9AP6. Hence differential RVs have "re-weighted" aperture data sets such that the color difference between aperture data sets disappears. This is also the reason why the periodogram difference between aperture data sets disappears if we include differential RVs in the model (see Figure 3). Figure 17. Lagrange stability of planets around τ Ceti. The shaded region represents the unstable orbits around a planet. The signals around 20, 49, 160, and 600 days are quantified simultaneously in the five-planet solution for the C1AP1 data set, and are denoted by red dots with error bars (measured at the 1% and 99% quantiles of posterior density). The shaded regions denote the parameter space of unstable orbits. Figure 19. MAP estimation of parameters varying with wavelength for the 9AP data sets. The error bars are determined by the 1% and 99% quantiles of the posterior densities. The MAP estimation of parameters for the 1AP1 data set and corresponding uncertainties are denoted by solid and dashed lines, respectively. The variation of d i (i 1 ,..., 8 Î { } ) with wavelength is step-like because two differential RVs with close wavelength ranges are derived from the same aperture data set but with opposite signs.