A publishing partnership

The following article is Open access

Color Difference Makes a Difference: Four Planet Candidates around τ Ceti

F. Feng, M. Tuomi, H. R. A. Jones, J. Barnes, G. Anglada-Escudé, S. S. Vogt, and R. P. Butler

Published 2017 September 5 • © 2017. The American Astronomical Society. All rights reserved.
The Astronomical Journal, Volume 154, Number 4 Citation F. Feng et al 2017 AJ 154 135 DOI 10.3847/1538-3881/aa83b4

Download Article PDF

DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

Article metrics

4307 Total downloads
Video abstract views

Dates

Received 2016 October 19
Revised 2017 July 27
Accepted 2017 July 31
Published 2017 September 5

Keywords

methods: numerical; methods: statistical; stars: individual (HD 10700); techniques: radial velocities

Journal RSS

Create or edit your corridor alerts

What are corridors?

Abstract

The removal of noise typically correlated in time and wavelength is one of the main challenges for using the radial-velocity (RV) method to detect Earth analogues. We analyze τ Ceti RV data and find robust evidence for wavelength-dependent noise. We find that this noise can be modeled by a combination of moving average models and the so-called "differential radial velocities." We apply this noise model to various RV data sets for τ Ceti, and find four periodic signals at 20.0, 49.3, 160, and 642 days, which we interpret as planets. We identify two new signals with orbital periods of 20.0 and 49.3 days while the other two previously suspected signals around 160 and 600 days are quantified to a higher precision. The 20.0 days 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 High Accuracy Radial Velocity Planet Searcher (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 RV technique to detect relatively weak signals.

Export citation and abstract BibTeX RIS

Previous article in issue

Next article in issue

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. 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⁻¹ 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⁻¹ 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 (Tuomi et al. 2013; 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⁻¹ (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 wavelength-dependent 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.

2. Radial-velocity Data of τ Ceti

In the European Southern Observatory archive, there are more than 9000 publicly available RVs measured by HARPS from 2003 June to 2013 September for HD 10700 as part of the observing programs 60.A-9036, Mayor, Comm, 072.C-0488, 072.C-0513, 074.D-0380, 075.C-0234, 075.D-0760, 076.C-0073, 077.C-0530, 078.C-0751, 078.C-0833, 079.C-0681, 081.C-0034, 082.C-0315, 083.C-1001, 084.C-0229, 085.C-0318, 086.C-0230, 087.C-0990, 088.C-0011, 089.C-0050, 090.C-0849, and 091.C-0936.

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,..., 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⁻¹), 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 wavelength-dependent 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.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** RV data and activity indicators of 1AP1 or HARPS (left panels) and Keck (right panels). The red lines denote the zero RV, and the blue line shows the epoch of JD 2453500, which is used to define the C1AP1 data set. The activity indices are normalized to the zero mean and unit standard deviation. The measurement uncertainties of the RV data are shown with error bars. There are 8880 and 752 data points in the 1AP1 and Keck data sets, respectively.
Download figure:
Standard image High-resolution image

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Data reduction and modeling process.
Download figure:
Standard image High-resolution image

3. 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 ${ \mathcal D }$ is the data, $P({M}_{i}| { \mathcal D })$ is the posterior, and $P({ \mathcal D }| {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.⁶ This is equivalent to the criterion, ${\rm{\Delta }}\mathrm{BIC}\gt 10$ , 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 ${ \mathcal D }$ , we need to estimate its parameters $\theta$ by calculating their posterior densities according to

where $P({ \mathcal D }| {\boldsymbol{\theta }},M)$ and $P({\boldsymbol{\theta }}| M)$ 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}_{k}| { \mathcal D },\,M)$ converges to a stationary distribution.

4. Differential Radial Velocities

4.1. 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 ${\rm{\Psi }}(t,\lambda )$ , the differential RV is defined by the difference between two aperture data sets,

where $v({t}_{i},{\lambda }_{j})$ is the RV at time t_i and wavelength ${\lambda }_{j}$ . If the wavelength difference is small enough, $D({t}_{i},{\lambda }_{j})$ is approximately proportional to the first partial derivative of Ψ with respect to ${\lambda }_{j}$ ,

Similarly, if the wavelength difference is a constant ( ${\rm{\Delta }}\lambda$ ) for all differential RVs, the second partial derivative of Ψ is approximately

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}_{\lambda }-1$ is the number of independent differential RVs, d_m characterizes the linear dependence of RV noise on differential RV at wavelength ${\lambda }_{m}$ , which is averaged over the differential RV, a is the acceleration caused by activity cycles or other wavelength-dependent noise⁷ 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

where ${f}_{k}({t}_{i})$ is the RV variation caused by the $k\mathrm{th}$ planet, and $\hat{{\rm{\Psi }}}$ is the estimation of the noise component. In the Keplerian component, K_k, ${\omega }_{k}$ , ${\nu }_{k}$ , and e_k are the amplitude, longitude of periastron, true anomaly, and eccentricity for the $k\mathrm{th}$ 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 ${\epsilon }_{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 $q\mathrm{th}$ -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 ${\lambda }_{j}$ 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

where ${\sigma }_{i,j}$ is the measurement noise at time t_i in the jth aperture data set, which has an averaged wavelength of ${\lambda }_{j}$ , ${s}_{J}({\lambda }_{j})$ 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.

Table 1. The Prior Distributions of Model Parameters

Parameter	Unit	Prior distribution	Minimum	Maximum
Each Keplerian Signal
K_j	m s⁻¹	$1/({K}_{\max }-{K}_{\min })$	0	$2\| v-\bar{v}{\| }_{\max }$
P_j	day	${P}_{j}^{-1}/\mathrm{log}({P}_{\max }/{P}_{\min })$	1	${t}_{\max }-{t}_{\min }$
e_j	⋯	${ \mathcal N }(0,0.1)$	0	1
${\omega }_{j}$	rad	$1/(2\pi )$	0	$2\pi$
M_0j	rad	$1/(2\pi )$	0	$2\pi$

Linear Trend and Jitter
a	m s⁻¹ yr⁻¹	$1/({a}_{\max }-{a}_{\min })$	$-365.24{K}_{\max }/{P}_{\max }$	$365.24{K}_{\max }/{P}_{\max }$
b	m s⁻¹	$1/({b}_{\max }-{b}_{\min })$	$-{K}_{\max }$	${K}_{\max }$
s_J	m s⁻¹	$1/({s}_{J\max }-{s}_{J\min })$	0	${K}_{\max }$

Moving Average
w	⋯	$1/({w}_{\max }-{w}_{\min })$	−1	1
τ	day	${P}_{j}^{-1}/\mathrm{log}({\tau }_{\max }/{\tau }_{\min })$	$1/({t}_{\max }-{t}_{\min })$	1

Activity Indices and Differential RVs
c_k	m s⁻¹	$1/({c}_{k\max }-{c}_{k\min })$	$-{c}_{k\max }$	${K}_{\max }/({I}_{X\max }-{I}_{X\min })$
d_m( $m\in \{1,\,\ldots ,\,{N}_{\lambda }-1\}$ )	m s⁻¹	$1/({d}_{m\max }-{d}_{m\min })$	$-{d}_{m\max }$	${K}_{\max }/({D}_{m\max }-{D}_{m\min })$

Note. The unit of c_k and d_m is m s⁻¹ 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-\bar{v}{| }_{\max }$ . The parameter characterizing the dependence of RV on activity indices is c_k, where k denotes the names of various indices.

Download table as: ASCII Typeset image

Nearly all RV data analysis has only considered the averaged data or 1AP1 is analyzed without accounting for wavelength-dependent 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.

4.2. 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⁸ ; Feng et al. 2017a; also see Baluev 2008 and Sü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 periodograms like the ones introduced by Baluev (2013, 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 wavelength-dependent 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.

Figure 3. Refer to the following caption and surrounding text. — **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.
Download figure:
Standard image High-resolution image

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.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Similar to Figure 3, but for the binned1AP1 (top left) and 1AP1 (bottom left) data sets and for their residuals (right panels) calculated by subtracting the 0P+MA(1)+2D prediction from the data.
Download figure:
Standard image High-resolution image

4.3. 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⁻¹ or deviate from the mean more than $5\sigma$ 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 $3\sigma$ . 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⁻¹, 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⁻¹ 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⁻¹. 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.

Figure 5. Refer to the following caption and surrounding text. — **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⁻¹, characterizing the level of instrumental bias. The mode, mean, standard deviation, and higher order moments of the posterior density are also shown.
Download figure:
Standard image High-resolution image

4.4. 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.

Table 2. BIC-Estimated BFs of Noise Models (Without Keplerian Component) for the C1AP1 Data Set

	MA(0)	MA(1)	MA(2)	MA(3)	MA(4)	MA(5)	MA(6)
0D	0.00	2234.55	2506.11	2606.66	2648.22	2649.77	2645.32
2D	411.11	2576.22	2840.77	2961.32	3013.88	3016.43	3011.98
5D	504.77	2596.88	2868.43	2989.98	3041.54	3041.09	3039.65
8D	537.43	2609.54	2883.09	3005.65	3053.20	3056.75	3053.31
17D	548.41	2596.52	2867.08	2994.63	3043.18	3044.74	3041.29

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.

Download table as: ASCII Typeset image

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.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Comparison between the periodograms for the residuals of the 1AP1 data set after subtracting the 0D+MA(4)+8D (left) and 0D+MA(4)+71D (right) prediction. The horizontal lines denote the FAPs of 0.1,1, and 0.001.
Download figure:
Standard image High-resolution image

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 wavelength-dependent 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.

4.5. 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 Figure 7. We see a steady increase of RV by 5 m s⁻¹ and a rapid decrease by 10 m s⁻¹ around the above epochs. These variations are huge compared with the sub-meter semi-amplitudes 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.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** RVs of all observation epochs (top) and for epochs around JD 2453283 (bottom left) and JD 2456195 (bottom right) for the 1AP1 data set. The bottom two panels show the RVs measured in the time spans denoted by the green lines in the top panel. The black and red dots with error bars represent binned RVs and binned residuals after subtracting the best-fitting noise model, respectively.
Download figure:
Standard image High-resolution image

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 intra-night 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.

5. 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.⁹ From now on, we only use the C1AP1 data set to identify signals and use 1AP1 and CCF to test the consistency of signals.

5.1. 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.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Logarithmic BF distributions of the samples drawn by the cold chains for the Keplerian (left) and circular (right) solutions of the 1P+MA(4)+8D model for the C1AP1 data set. There are 10,000 samples randomly drawn from the each cold chain. The periods of some local maxima are denoted by red lines and numbers. If two maxima are close to each other, we only show the one with higher posterior. The BIC-estimated logarithmic BF threshold of 5 is denoted by the dashed line in each panel.
Download figure:
Standard image High-resolution image

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.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Logarithmic likelihood ratio of the 1P+MA(4)+8D model and the 0P+MA(4)+8D model based on the 3.8 million posterior samples drawn by hot chains for the period intervals of days (left) and days (right) for the C1AP1 data set. Each period range is divided into 1000 bins, and the maximum posterior is determined for each bin. The posterior distributions shown here are linear interpolations of these posteriors. The 10%, 1%, and 0.1% quantiles of the maximum logarithmic likelihood ratio are denoted by dashed lines. Note that the likelihood ratio determined by hot chains does not represent the real one.
Download figure:
Standard image High-resolution image

$[13.5,20.9]$ — **Figure 9.** Logarithmic likelihood ratio of the 1P+MA(4)+8D model and the 0P+MA(4)+8D model based on the 3.8 million posterior samples drawn by hot chains for the period intervals of days (left) and days (right) for the C1AP1 data set. Each period range is divided into 1000 bins, and the maximum posterior is determined for each bin. The posterior distributions shown here are linear interpolations of these posteriors. The 10%, 1%, and 0.1% quantiles of the maximum logarithmic likelihood ratio are denoted by dashed lines. Note that the likelihood ratio determined by hot chains does not represent the real one.
Download figure:
Standard image High-resolution image

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.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Similar to Figure 8, but for cold chains for the C1AP1 and CCF data sets.
Download figure:
Standard image High-resolution image

5.2. 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 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.

Table 3. MAP Estimation of the Parameters of Signals for the CCF and C1AP1 Data Sets and for Both Keplerian (Left Columns) and Circular Solutions (Right Columns)

	Keplerian		Circular
	CCF	C1AP1	CCF	C1AP1
P₁	1007.66	995.78	1015.25	942.68
	597.02	636.13	598.55	619.39

K₁	0.70	0.39	0.30	0.26
	0.68	0.35	0.43	0.36

e₁	0.46	0.30	⋯	⋯
	0.58	0.16	⋯	⋯

P₂	164.38	161.40	164.28	162.99
	164.22	162.87	164.72	163.42

K₂	0.53	4.65	0.44	0.49
	0.48	0.55	0.47	0.51

e₂	0.26	0.96	⋯	⋯
	0.18	0.18	⋯	⋯

P₃	20.03	20.03	20.03	20.01
	20.04	20.00	20.03	20.03

K₃	0.47	0.42	0.41	0.46
	0.49	0.49	0.39	0.43

e₃	0.18	0.11	⋯	⋯
	0.25	0.06	⋯	⋯

P₄	49.47	49.29	49.36	49.37
	49.38	49.41	49.48	49.50

K₄	0.45	0.48	0.38	0.32
	0.42	0.39	0.42	0.42

e₄	0.28	0.29	⋯	⋯
	0.22	0.23	⋯	⋯

P₅	91.59	102.43	92.29	91.97
	91.97	91.34	91.65	91.72

K₅	0.52	0.46	0.27	0.28
	0.32	0.45	0.32	0.30

e₅	0.53	0.15	⋯	⋯
	0.27	0.40	⋯	⋯

P₆	⋯	⋯	6.63	⋯
	⋯	⋯	6.63	31.68

K₆	⋯	⋯	0.28	⋯
	⋯	⋯	0.25	0.30

e₆	⋯	⋯	⋯	⋯
	⋯	⋯	⋯	⋯

P₇	⋯	⋯	14.21	⋯
	⋯	⋯	14.22	⋯

K₇	⋯	⋯	0.27	⋯
	⋯	⋯	0.19	⋯

e₇	⋯	⋯	⋯	⋯
	⋯	⋯	⋯	⋯

Note. The parameters shown are the period (P in days), semi-amplitude (K in m s⁻¹), and eccentricity (e) of the signals. For circular solutions, only semi-amplitude and period are shown. For each signal, the parameters are estimated in combination with the ∼1000 day and with the ∼600 day signals are shown in upper and lower entry, respectively, for each parameter and data set. The solutions are obtained based on at least 4 million MCMC samples.

Download table as: ASCII Typeset image

Table 4. Logarithm BFs Estimated by the BIC for Models with Various Numbers of Planets with Respect to the Noise Model

Number of Planets	1	2	3	4	5	6
Keplerian solution (∼1000 days)	46.6	89.7	191	219	⋯	⋯
Keplerian solution (∼600 days)	63.7	85.3	126	152	⋯	⋯
Circular solution (∼1000 days)	11.6	59.1	95.3	203	112	⋯
Circular solution (∼600 days)	26.8	72.6	104	118	125	134

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.

Download table as: ASCII Typeset image

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 high-eccentricity 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.

Figure 11. Refer to the following caption and surrounding text. — **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).
Download figure:
Standard image High-resolution image

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 long-period 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.

5.3. 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⁻¹ 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⁻¹. 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⁻¹ 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 parameters varying with chunks. This suggests that the signals we have identified are probably not caused by activity cycles.

5.4. 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,

where

where n is a natural number. The zero moment is simply the sum of flux values, ${M}_{0}={\sum }_{i=1}^{N}{f}_{i}(t)$ . The central moments are defined by

To minimize excess instrumental noise, we use the central moments together with M₀ because the non-central moments contain a certain amount of information from Keplerian signals. Although the central moments attempt to minimize the mean velocity contribution by subtracting M₁, M₁ 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₀, M_3c, M_4c, and M_5c,¹⁰ 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₀, M_3c, M4c, and M_5c are around −0.17 ± 0.02 m s⁻¹, −0.62 ± 0.02 m s⁻¹, 0.16 ± 0.02 m s⁻¹, and −0.28 ± 0.02 m s⁻¹, respectively.¹¹ 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₁ 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.

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Correlation between RVs and M₀, M_3c, M_4c,and M_5c for the C1AP1 data set. The Pearson correlation coefficient is shown in the top right corner.
Download figure:
Standard image High-resolution image

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 longer than one epoch. The Pearson correlation coefficients between RVs and proxies for the epochs around JD 2453283 are −0.40, −0.36, 0.28, −0.02, and 0.02 for the M₀, M_c3, M_c4, Mc5, and cC3AP2-1, respectively, while the coefficients are 0.06, −0.73, 0.43, −0.50, and 0.32 for the other epochs. 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.

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Normalized binned RVs, moments, and c3AP2-1 data set (CDS) around epochs of JD 2453283 and JD 2456195 for the 1AP1 data. The binning windows for the left and right panels are 0.05 days and 0.5 days, respectively. The c3AP3-2 data set around JD 2453283 is set to the value of nearby epochs since no calibration stars are available around this epoch. The mean of each time series is shown by the horizontal dashed line while all time series are shifted up/down for visualization.
Download figure:
Standard image High-resolution image

Although the moments correlate with RVs up to 0.5 m s⁻¹, 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 signal-to-noise models.

6. Planetary Candidates

6.1. 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⁻¹,¹² which is beyond the detection ability of Keck.

Figure 14. Refer to the following caption and surrounding text. — **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.
Download figure:
Standard image High-resolution image

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, ${\theta }_{1 \% }$ and ${\theta }_{99 \% }$ , are determined at the cumulative marginalized posterior probability of 1% and 99%. Specifically, they are determined by ${\int }_{{\theta }_{\min }}^{{\theta }_{1 \% }}P(\theta | { \mathcal D })d\theta =1 \%$ and ${\int }_{{\theta }_{\min }}^{{\theta }_{99 \% }}P(\theta | { \mathcal D })d\theta =99 \%$ , where ${\theta }_{\min }$ is the minimum value of parameter θ.¹³ 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}_{\oplus }$ . 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⁻¹, which is close to the 0.1 m s⁻¹ limit required for detecting Earth analogues around solar analogues.

Table 5. The MAP Estimation of the Parameters for Four Signals Detected in the C1AP1 Data Set

Parameters	τ Ceti g	τ Ceti h	τ Ceti e	τ Ceti f
P (days)	20.00 [19.99, 20.02]	49.41 [49.31, 49.49]	162.87 [162.41, 163.95]	636.13 [588.44, 647.83]
K (m s⁻¹)	0.49 [0.38, 0.56]	0.39 [0.33, 0.54]	0.55 [0.46, 0.68]	0.35 [0.23, 0.45]
e	0.06 [0.00, 0.19]	0.23 [0.08, 0.39]	0.18 [0.04, 0.36]	0.16 [0.00, 0.23]
ω(rad)	6.90 [5.42, 7.53]	0.13 [−0.71, 0.78]	0.39 [−0.41, 1.68]	2.09 [0.86, 2.81]
M₀(rad)	7.04 [6.67, 8.65]	−1.27 [−1.84, −0.52]	6.21 [5.27, 7.37]	−0.68 [−1.32, 0.30]
$m\sin i$ ( ${M}_{\oplus }$ )	1.75 [1.35, 2.00]	1.83 [1.57, 2.51]	3.93 [3.29, 4.76]	3.93 [2.56, 4.98]
a(au)	0.133 [0.131, 0.134]	0.243 [0.240, 0.246]	0.538 [0.532, 0.544]	1.334 [1.290, 1.351]

Note. The uncertainties of parameters are denoted by the lower and upper bounds of the intervals, which are determined at the 1% and 99% quantiles of the posterior density. To estimate the semimajor axis a and the minimum mass of planet $m\sin i$ , we set the mass of τ Ceti to be $0.783\pm 0.012{M}_{\odot }$ (Teixeira et al. 2009). To be consistent with the names used by MT13, the new signals at periods of 20 and 49 days are named τ Ceti g and h.

Download table as: ASCII Typeset image

If we regard the measurements in a 15 minute time bin as independent observations (e.g., Mayor et al. 2003 and O'Toole et al. 2008), there would be 662 observations for the C1AP1 data set. The $K/N=K/{{RV}}_{\mathrm{rms}}\times \sqrt{{N}_{\mathrm{obs}}}$ ratio for signals with K = 0.3 m s⁻¹ is 7.3, where ${{RV}}_{\mathrm{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.

Figure 15. Refer to the following caption and surrounding text. — **Figure 15.** The phase-folded data and the best model prediction for each signal for the C1AP1 data set. The data and the binned version are shown as small and big dots. The MAP estimation of the four signals are based on 8 million posterior samples drawn by a cold chain.
Download figure:
Standard image High-resolution image

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.

Figure 16. Refer to the following caption and surrounding text. — **Figure 16.** 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.
Download figure:
Standard image High-resolution image

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.

6.2. 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⁻¹. These values are lower than the values reported by MT13, which are from 0.58 to 0.75 m s⁻¹. 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 wavelength-dependent 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.

6.3. 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.

Figure 17. Refer to the following caption and surrounding text. — **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.
Download figure:
Standard image High-resolution image

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 (MacGregor et al. 2016). The dynamical simulations of the planets reported by MT13 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}_{\odot }$ , 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.

7. 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 wavelength-dependent 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⁻¹ 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 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/.

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⁻¹. 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⁻¹. 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.

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⁻¹ 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).

Figure 18. Refer to the following caption and surrounding text. — **Figure 18.** Correlation between noise model parameters based on 1000 samples randomly drawn from 4 million posterior samples of the 4P+MA(4)+8D model for the C1AP1 data set. In this Keplerian solution, the signals have periods of 49.5, 91.7, 20, and 597 days. The Pearson correlation coefficient is shown in the top left corner for each panel. The parameters of c_s, c_b, c_f, and c_c are the linear coefficients of S-index, BIS, FWHM, and cC3AP2-1. The other parameters are linear coefficients of differential RVs derived from 9AP aperture data sets. All parameters are in units of m s⁻¹.
Download figure:
Standard image High-resolution image

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 19. Refer to the following caption and surrounding text. — **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\in \{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.
Download figure:
Standard image High-resolution image

**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\in \{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.
Download figure:
Standard image High-resolution image

Footnotes

6
Although the comparison of BF estimators performed by Feng et al. (2016) is for M dwarfs, the BIC is found to be rather conservative, and thus is appropriate for τ Ceti, which is a quiet star.
7
For example, the atmospheric absorption and/or scattering of the star light could be wavelength dependent, leading to wavelength-dependent noise in RVs. Considering that the linear trend in the model may be caused by stellar activity and instrumental noise, we put it in the noise component of the RV model.
8
The code for calculating GLST is available at https://github.com/phillippro/agatha and a corresponding online app is at http://www.agatha.herts.ac.uk.
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.
10
We do not use M_2c because it is equivalent to FWHM.
11
These values are the mean and standard deviation of the posterior samples drawn by cold chains, and thus are different from the Pearson coefficients shown in Figure 12.
12
The mean semi-amplitudes of the Keplerian signals are equal to the semi-amplitudes determined in the circular solutions, which are shown in the right columns of Table 3.
13
Actually, the formula expressed here should be expressed in a discrete form because we calculate the intervals using MCMC samples.

Please wait… references are loading.

Color Difference Makes a Difference: Four Planet Candidates around τ Ceti

Article metrics

Share this article

Dates

Abstract

1. Introduction

2. Radial-velocity Data of τ Ceti

3. Data Analysis Methods

4. Differential Radial Velocities

4.1. Wavelength-dependent Noise Models

4.2. Removing Wavelength-dependent Noise

4.3. Modeling Instrumental Noise

4.4. Choosing the Optimal Noise Model

4.5. Performance of the Goldilocks Noise Model

5. Keplerian Signals

5.1. Primary Signals

5.2. Comparing Signals Detected in Different Data Sets

5.3. Curse of Eccentricity

5.4. Instrumental Bias of HARPS

6. Planetary Candidates

6.1. Four-planet System

6.2. Comparison with Previous Studies

6.3. Dynamical Stability and Habitability Analysis

7. Discussions and Conclusions

Appendix A: Correlation between Noise Proxies

Appendix B: Wavelength dependence of Model Parameters

Footnotes