THE HARPS-TERRA PROJECT. I. DESCRIPTION OF THE ALGORITHMS, PERFORMANCE, AND NEW MEASUREMENTS ON A FEW REMARKABLE STARS OBSERVED BY HARPS

Guillem Anglada-Escudé; R. Paul Butler

doi:10.1088/0067-0049/200/2/15

1. INTRODUCTION

The Doppler technique has been the most successful method of detecting and confirming the presence of extrasolar planets around nearby stars (Mayor & Queloz 1995; Marcy & Butler 1996). The stability of the spectrographs and the data analysis techniques used to obtain precision radial velocity (RV) measurements have been steadily improving during the last 16 years of exoplanet discoveries. In particular, the vast majority of candidates detected via Doppler spectroscopy have been obtained using two approaches: the gas cell method and the stabilized spectrograph approach.

Because an absorption cell can be installed at low cost on any general-purpose echelle spectrometer, it is the most broadly used technique. In this approach, a glass cell filled with iodine gas at low pressure is inserted just before the entrance slit of the spectrometer. Because the stellar spectrum is imprinted with iodine prior to entering the spectrograph, the cell spectrum tracks the same instrumental distortions suffered by the stellar lines. Thanks to this, the wavelength solution and the instrumental profile can be fitted simultaneously to the Doppler shift of the star, allowing nominal RV precisions at the level of 1–2 m s⁻¹. The data analysis method required to extract such precise measurements consists of forward-modeling the spectrum of iodine multiplied by a high signal-to-noise ratio (S/N) template of the star and convolving the product with a parameterized model of the instrumental profile. This method was pioneered by Butler et al. (1996) and has been used to infer the presence of more than 300 exoplanets. Due to the large number of parameters involved, this technique is computationally very intensive. Also, it requires a high-S/N realization of the real stellar spectrum at higher resolution than the actual observations. Because obtaining a high-resolution perfectly calibrated spectrum of the star is usually not possible, obtaining a realistic template from the deconvolution of observations without iodine is a key element of the method and is thought to be one of its major limitations (e.g., see supplementary material in Howard et al. 2010).

The other leading method of obtaining precision RV measurements consists of building a fiber-fed, very stable spectrograph, which is then wavelength-calibrated using an external source such as a Th/Ar emission lamp. This technique has been developed and refined by the group lead by M. Mayor (hereafter "the Geneva group") over the past 20 years. This approach provided the first clear detection of an extrasolar planet around a solar-like star (Mayor & Queloz 1995) using the ELODIE spectrograph at the Observatoire de Haute-Provence in France. The interested reader can find the basic elements used in the construction of a fiber-fed stabilized spectrograph in Baranne et al. (1996). While the gas cell technique requires the simultaneous adjustment of the wavelength scale, instrumental profile, and Doppler offset of the star, the stabilized spectrograph approach allows each problem to be tackled separately. First, the instrumental profile is constant by design thanks to the use of an image scrambling system coupled with optical fibers. In the design by Baranne et al. (1996), the wavelength calibration is obtained using two optical fibers closely packed together at the entrance of the spectrometer. Therefore, they follow almost identical optical paths within the spectrometer optics tracking almost the same optical distortions. The first fiber (the science fiber) is illuminated with a wavelength calibration source (e.g., Th/Ar lamp) at the beginning of the night, and the nominal absolute wavelength solution relative to that source is obtained (Pepe et al. 2002). The second fiber is fed with the same calibration source during the calibration of the science fiber and during the science observations. Although there might be a significant RV offset between fibers, they are similarly affected by changes in the instrument and, as a consequence, they share similar intra-night wavelength drifts. This second fiber is only used to monitor intra-night changes in the instrument, and the information it provides is an RV drift to be added to the measured RVs. Finally, the RV measurement on a fully calibrated spectrum is obtained by the so-called cross-correlation function method (CCF). CCF is based on multiplying the observed spectrum by a weighted binary mask. The binary mask is different from zero on the nominal positions of the stellar lines, and each non-zero chunk is weighted according to the relative depth of the stellar line against the local continuum. This binary mask is then Doppler-shifted and the CCF is evaluated again. The minimum of the CCF as a function of the Doppler offset is the desired RV measurement. To achieve higher precision, the shape of the CCF is centroided using a Gaussian profile. The CCF method and some details on how the binary masks are obtained are outlined in Queloz (1995) and Pepe et al. (2002).

The flagship stabilized spectrometer built by the Geneva group is the High Accuracy Radial velocity Planet Searcher (HARPS) installed on the 3.6 m telescope at the European Southern Observatory (ESO) site in La Silla, Chile (Pepe et al. 2003). Thanks to a careful design and construction (a vacuum-sealed tank, high mechanical stability, and accurate temperature control), the wavelength solution and instrumental profile are very stable. For example, intra-night Doppler drifts as measured by the calibration fiber are typically smaller than 0.5 m s⁻¹. Long exposures (t >200 s) using the calibration fiber cause leaks of the Th/Ar spectrum on the science spectrum, so the calibration fiber is only used when observing very bright targets and when extreme RV precision is required. The list of planets detected by HARPS is long and varied as can be seen in the 34 papers of the series "The HARPS Search for Southern Extra-solar Planets." Instead of citing all of them, we refer the interested reader to the latest HARPS results presented in Pepe et al. (2011) and Mayor et al. (2011). HARPS has demonstrated an RV stability at the level of 1 m s⁻¹ on timescales of several years.

Despite these impressive results, it is known that the CCF method implemented in the HARPS Data Reduction Software (DRS) is suboptimal in the sense that it does not exploit the full Doppler information on the stellar spectrum (e.g., see Queloz 1995; Pepe et al. 2002). We asked ourselves if a least-squares approach where the observed stellar spectrum is matched to a high-S/N template could be used to extract higher RV precision on HARPS observations. Given that a number of stabilized spectrographs are under construction (e.g., HARPS-North, ESPRESSO/ESO, Carmenes/CaHa) and that there is significant investment in hardware development to achieve higher RV stability, it is important that the used data analysis methods are as optimal as possible. Also, precision RVs are difficult to reproduce, and, given that the stakes are ambitious (detection of potentially habitable worlds), the use of different RV measurement improves the detection confidence of low-amplitude signals. In this work we derive from first principles the algorithms of the template-matching technique to be applied on stabilized spectrographs and implement it to public HARPS observations. HARPS is an ESO instrument and, as such, all the data obtained from it become publicly available after a proprietary period of a few months (or years). Since 2011 January, the data products derived from the HARPS-ESO DRS (wavelength-calibrated spectra, but also CCF-measured RVs among others) are publicly available through a dedicated Web page in the ESO Web site.³ All the HARPS data used in this work have been obtained from there.

The core of our project is our software tool called HARPS-TERRA, where TERRA stands for Template-Enhanced Radial velocity Re-analysis Application. HARPS-TERRA handles the full process of unpacking the HARPS-ESO archive files, generation of a high-S/N template, and obtaining the final RV measurement in a single command line call. This software is custom made, fully coded in Java, and can run on any machine supporting a Java Run-time Environment 1.6 or newer.

In Section 2, we derive and describe the basic algorithms used to obtain RV measurements from HARPS-reduced spectra. In Section 3, we investigate the optimal RV extraction parameters and define the standard setup for obtaining optimal RV precision with minimum human intervention on G, K, and M dwarfs. Section 4 shows the performance of HARPS-TERRA on a few representative data sets, demonstrating a significant improvement in precision, especially with M dwarfs. In the same section, we use new RV measurements obtained with HARPS-TERRA to discuss planetary systems proposed on a few representative stars observed by HARPS.

HARPS-TERRA is still in development but can be distributed for particular applications upon request. Given that other HARPS programs could benefit from using it (e.g., asteroseismology programs, binary stars), we plan a public release of the tool in the near future.

2. DESCRIPTION OF THE ALGORITHMS

The proposed algorithm is based on minimizing the differences of the observed spectrum against a parameterized template. In our particular implementation, we first use the higher-S/N observation as a preliminary template. In a second iteration, a very high S/N template is obtained by co-adding all the observations and the RVs are measured again. Note that the template will be already convolved with the instrumental profile so this method (like the CCF technique) relies on the long-term stability of the instrumental profile and wavelength calibration strategy.

HARPS is a cross-dispersed echelle spectrograph. Therefore, the stellar spectrum is split in diffraction orders over the detector (also called echelle apertures or apertures). Each aperture has to be extracted very carefully in a complex process from the raw CCD images. The wavelength calibration is usually made on a nightly basis using a standard set of calibration frames (i.e., flat, darks, Th/Ar lamps) that are taken once at the beginning of the night. Thankfully, all this extraction and the corresponding wavelength solution are efficiently implemented by HARPS DRS⁴ and will not be discussed here. The wavelength-calibrated spectra are provided in conveniently formatted "fits" files described in the HARPS DRS manual. As of 2011 November, HARPS-TERRA is designed to work with the output files of HARPS DRS v3.5 but should be easily adapted to future updates of HARPS DRS and/or other stabilized spectrographs. An overview of the HARPS spectral format is given in Table 1. As described in the Introduction, the secondary calibration fiber is typically only used on bright stars to achieve maximum precision. If the observer used this secondary fiber, such a drift is also provided by HARPS DRS and will be added to the final RV measurement of each echelle aperture. In what follows we treat each echelle aperture as an independent spectrum and the final RV measurement will be a weighted mean of the RVs measured across all such apertures.

Table 1. HARPS Spectral Parameters Relevant to This Work

Parameter	Value
Spectral resolution λ/δλ at 5500 Å	120,000
Number of diffraction orders (or echelle apertures)	72
Pixels in each aperture	4096
Wavelength range (Blue CCD)	3780–5300 Å
Wavelength range (Red CCD)	5330–6910 Å
Sampling per resolution element	4.1 pixel

Download table as: ASCII Typeset image

Let us assume that we have a very high S/N, wavelength-calibrated spectrum of the star (template). Given a wavelength-calibrated observation and a stable instrumental profile, the observed spectrum differs from the template only by a Doppler shift (due to Earth's motion around the Sun and/or the presence of companions) and a flux normalization function across each aperture. We observed that this flux normalization is time-dependent at the few percent level most likely due to observational and instrumental effects such as atmospheric differential refraction, differential absorption, or telescope tracking errors. Let us define the difference R between the template F and the observed flux f at each wavelength λ as

$\begin{eqnarray} &&R\left[\lambda ;\hat{\boldsymbol{\alpha }}\right] = F\left[ \alpha _v \lambda \right]- f\left[\lambda \right] \sum _{m=0}^{M} \alpha _m \left(\lambda -\lambda _c\right)^{m}, \end{eqnarray} \tag{ 1 }$

where λ is the wavelength of the observation transformed to the solar system barycenter reference frame. The set of free parameters is represented by $\hat{\boldsymbol{\alpha }}= [\alpha _v, \alpha _0, \ldots , \alpha _M]$ . The first term on the right side is the template evaluated at α_vλ, where α_v is the Doppler factor in which we are interested. To derive a differential RV measurement from this Doppler factor, we use the very simple expression α_v = 1 − v_r/c (Stumpff 1979), where c is the speed of light and v_r is the relative RV between the observer and the star in which we are interested. The second term is the observed flux f multiplied by an M-degree polynomial accounting for the flux normalization across each aperture (also called blaze function correction), and λ_c is the central wavelength of a given aperture. Note that we could apply this flux normalization to the template F instead of to the observed flux f. However, this would couple the flux normalization coefficients α_m to the Doppler factor α_v in a nonlinear fashion, which, from a numerical point of view, is an undesirable complication.

Now this difference R can be Taylor-expanded around some nominal values for the parameters $\hat{\boldsymbol{\alpha }}_{(0)}$ in powers of the parameter increments $\delta \hat{\boldsymbol{\alpha }}$ as

$\begin{eqnarray} R[\lambda ;\hat{\boldsymbol{\alpha }}] &\simeq & R[\lambda ;\hat{\boldsymbol{\alpha }}_{(0)}] + \left. \frac{\partial R}{\partial \alpha _v} \right|_{\lambda ;\,\hat{\boldsymbol{\alpha }}_{(0)}} \delta \alpha _v + \sum _{m=0}^{M} \left. \frac{\partial R}{\partial \alpha _m} \right|_{\lambda ;\,\hat{\boldsymbol{\alpha }}_{(0)}} \delta \alpha _m \nonumber\\ \end{eqnarray} \tag{ 2 }$

$\begin{eqnarray} \left. \frac{\partial R}{\partial \alpha _v} \right|_{\lambda ;\,\hat{\boldsymbol{\alpha }}_{(0)}} &=& \left. \lambda \frac{d F}{d (\alpha _v\lambda)} \right|_{\alpha _{v(0)}\lambda } \end{eqnarray} \tag{ 3 }$

$\begin{eqnarray} \left. \frac{\partial R}{\partial \alpha _m} \right|_{\lambda ;\,\hat{\boldsymbol{\alpha }}_{(0)}} &=& -f[\lambda]\left(\lambda -\lambda _c\right)^m \,. \end{eqnarray} \tag{ 4 }$

The weighted sum of these R over all the observed wavelengths λ_i (or pixels) is defined as

$\begin{eqnarray} &&\chi ^2 = \sum ^{N_{{\rm pixel}}}_{i=1} \omega _i R[\lambda _i;\hat{\boldsymbol{\alpha }}]^2 \end{eqnarray} \tag{ 5 }$

and is the quantity to be minimized. The partial derivative of χ² with respect to each increment δα_v, δα₀, δα₁, ...δα_M equated to 0 generates the system of equations for such increments (the so-called normal equations), which can then be solved using standard matrix techniques. The quantities ω_i are the weights assigned to each λ_i, and their precise value will be discussed later. The resulting system of equations reads

$\begin{eqnarray} &&\sum _{i}^{N_{{\rm pixel}}} \left. \omega _i \frac{\partial R}{\partial \alpha _v} \frac{\partial R}{\partial \alpha _v} \right|_{\lambda _i;\,\hat{\boldsymbol{\alpha }}_{(0)}} \, \delta \alpha _v + \cdots + \sum _{i}^{N_{{\rm pixel}}} \left. \omega _i \frac{\partial R}{\partial \alpha _v} \frac{\partial R}{\partial \alpha _M} \right|_{\lambda _i;\,\hat{\boldsymbol{\alpha }}_{(0)}} \, \nonumber\\ &&\times\delta\alpha _M= -\sum _{i}^{N_{{\rm pixel}}} \left. \omega _i R \frac{\partial R}{\partial \alpha _v} \right|_{\lambda _i;\,\hat{\boldsymbol{\alpha }}_{(0)}} \end{eqnarray} \tag{ 6 }$

$\begin{eqnarray} &&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ldots \nonumber \\ &&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ldots \nonumber \\ &&\sum _{i}^{N_{{\rm pixel}}} \left. \omega _i \frac{\partial R}{\partial \alpha _M} \frac{\partial R}{\partial \alpha _v} \right|_{\lambda _i;\,\hat{\boldsymbol{\alpha }}_{(0)}} \, \delta \alpha _v + \cdots + \sum _{i}^{N_{{\rm pixel}}} \left. \omega _i \frac{\partial R}{\partial \alpha _M} \frac{\partial R}{\partial \alpha _M} \right|_{\lambda _i;\,\hat{\boldsymbol{\alpha }}_{(0)}} \nonumber\\ &&\times \delta \alpha _M= -\sum _{i}^{N_{{\rm pixel}}} \left. \omega _i R \frac{\partial R}{\partial \alpha _M} \right|_{\lambda _i;\,\hat{\boldsymbol{\alpha }}_{(0)}}, \end{eqnarray} \tag{ 7 }$

where the values of R are obtained using Equation (1) and the partial derivatives are computed using Equations (3) and (4). Because the Doppler factor α_v is nonlinear, this system of equations has to be solved iteratively. Each iteration consists of (1) computing the normal equations using the current values of the parameters $\hat{\boldsymbol{\alpha }}_{(0)}$ , (2) solving for the parameter increments $\delta \hat{\boldsymbol{\alpha }}$ , and (3) updating the parameter values $\hat{\boldsymbol{\alpha }}_{{\rm new}} = \hat{\boldsymbol{\alpha }}_{0}+\delta \hat{\boldsymbol{\alpha }}$ . These equations can be written in matrix form as

$\begin{eqnarray} &&A_{lk}\delta \alpha _{k} = b_l, \end{eqnarray} \tag{ 8 }$

where A_lk is the M+1 × M+1 matrix of coefficients multiplying the parameter increments δα_k in the normal Equations (6) and (7). A_lk is equivalent to the Hessian matrix of the χ² and is sometimes called the curvature matrix. By straightforward error propagation (Press et al. 1992) one finds that, when the solution converges to the χ² minimum, the formal uncertainties in each free parameter can be derived from the inverse of the $\mathbf {A}$ matrix as

$\begin{eqnarray} &&\sigma _{\alpha _k} = \sqrt{(\mathbf {A}^{-1})_{kk}}\, , \end{eqnarray} \tag{ 9 }$

where $\mathbf {A}^{-1}$ is usually called the covariance matrix. Even though this prescription ignores correlation between parameters, it still provides a good estimate of the formal precision of the measurement if the weights ω_i are properly estimated. The formal uncertainty in the Doppler coefficient δα_v multiplied by the speed of light (σ_v = c δα_v) is used to compute the initial weight of each echelle aperture in the computation of the mean epoch RV. Therefore, one needs to be sure that these formal uncertainties are as realistic as possible (see Section 2.1.2).

The RV measurement obtained on each aperture, its formal uncertainty, and other information about the fit (flux normalization coefficients, number of iterations, number of masked pixels, rms of the fit, etc.) are stored in a file for further processing (see Section 2.2). At this point, no assumption is made about the nature of the star and all the apertures are processed irrespective of their S/N or overall quality.

2.1. Algorithm Implementation Details

Even though the basic algorithm has already been outlined, a number of technical details must be carefully addressed to reach maximal RV precision in an efficient numerical fashion. Here we describe those technical aspects that require special attention and provide the practical solutions we have adopted.

2.1.1. Template Interpolation

Until now, we have assumed a theoretical continuous and differentiable function for the template. In reality, we need to build the template from the same observations. Since the observed wavelengths do not necessarily coincide with the wavelengths where the template is sampled, F and its derivatives must be obtained through interpolation. In our implementation, the interpolation of the discretely sampled template is obtained using a cubic spline (Press et al. 1992). Cubic splines have the property of producing a continuous version of the derivative, which is critical for the convergence properties of the nonlinear least-squares algorithm. The details on how such a template is generated are given in Section 2.3.

2.1.2. Pixel Weighting

Each weight ω_i should be a function of the uncertainty in R_i. Assuming a very high S/N template and Poisson statistics, the uncertainty in the flux of each pixel is just $\sigma _i=\sqrt{f_i}$ . However, we found that assuming such uncertainty produced suboptimal results because (1) it tends to overweight pixels with high flux and little real Doppler information (continuum) and (2) Poisson statistics never apply in realistic situations (e.g., pixel nonlinearities, imperfect flat-fielding, unmasked telluric features, imperfect template). Instead, we re-scale the Poisson uncertainty of each pixel using an empirical approach. In the first iteration of Equations (6) and (7), Poisson statistics are assumed on each pixel so that $\omega _{i(0)}=\sqrt{f_i}$ . After that, we compute the rms of the R_i differences over all the pixels and call it σ_〈R〉. If Poisson statistics are applied, σ_〈R〉 should be equal to the square root of the mean pixel flux ( $\sqrt{\langle f\rangle}$ ). However, if there are additional sources of noise, $\sigma _{\langle R\rangle}/\sqrt{\langle f\rangle}= \kappa$ , where κ > 1. The weights to be used in the next iteration are therefore re-scaled with this quality factor κ as ω_i(new) = (κ²f_i)⁻¹. The κ value of the last least-squares iteration is also stored in a file and provides useful information to assess the level of systematic noise on each echelle order.

2.1.3. Telluric Masking and Outlier Filtering

The weights are also used to mask those pixels with undesirable properties. That is, any ω_i on wavelengths coincident with telluric absorption features deeper than 1% is set to 0. A synthetic spectrum of the atmosphere was used to identify such telluric lines and generate a list of wavelengths (or telluric mask) to be avoided. The same telluric mask is applied to all the observations.

Because of the barycentric motion of the Earth (∼ ± 15 km s⁻¹), the stellar spectrum at the border of each echelle order is not always present in all the epochs. This going in and out of the borders can generate small systematic RV offsets correlated with the Earth barycentric motion. Even though a more sophisticated solution could be adopted and for the sake of simplicity, the weights of pixels within 0.8 Å of each aperture border are also set to 0. Although some Doppler information is lost in the process, the flux at the borders is significantly lower than at the center of the apertures, and some extra random noise is preferred over correlated noise with a one-year periodicity. We tested less restrictive border masking and obtained almost identical results.

After the first least-squares iteration is obtained, the weights of possible outliers are also set to 0. A pixel is considered an outlier if it has a residual four times larger than the empirically determined σ_〈R〉 (i.e., 4σ clipping). This threshold is arbitrary, but we found that it does a good job removing pixel outliers (e.g., cosmic ray hits), while preserving most of the well-behaved ones.

2.1.4. Parameter Initialization and Barycentric Correction

The described algorithm is a nonlinear least-squares solver, and, as such, it requires an initial guess for the values of the free parameters. The initial value for the flux polynomial coefficient α₀ is set to the mean pixel flux 〈f〉, and all the other α_m are set to 0. Note that the real flux normalization function (or blaze function) of the template is unknown. However, we are only interested in correcting relative changes compared to the template, so a low-order polynomial correction should be sufficient given a stabilized instrument such as HARPS. Section 3.1 investigates the optimal degree required for this polynomial.

All the analysis described above is done in barycentric wavelengths. The correction from observed to barycentric wavelengths is also implemented using the aforementioned simplistic recipe for the Doppler factor using the projected RV of the observer as computed by Stumpff (1980) (also provided in the header of the HARPS DRS files). Therefore, as a first guess, α_v is initialized to 1. Before the least-squares adjustment begins, α_v is initialized by finding the approximate χ² minimum using Equation (5) in a bisection algorithm with α_v as the independent variable. The obtained RV from this bisection is typically within 10 m s⁻¹ of the final RVs, minimizing the required number of computationally expensive least-squares iterations.

Note that the simplistic definition of the Doppler factor we are using is only accurate to first order on v/c. Since we are only interested in differential RV measurements, one can show that the second-order terms (∼v²/c² but also gravitational terms proportional to G) are mostly constant and that this expression for α_v is good at a few tens of cm s⁻¹. While higher-order terms should certainly be taken into account, such refinement is useless unless a full revision of the barycentric correction model is implemented. Critical parts in this model that need to be revised are (1) the use of an up-to-date model for Earth rotation, (2) the use of the most recent version of the ephemeris for Earth (e.g., DE405 or newer; Standish 1998), and (3) inclusion of all the relativistic corrections to the first post-Newtonian order (gravitational and kinematic). Even though this model might improve some of the RVs presented here, further work is required to produce a functional and reliable code. A tool to provide barycentric Doppler corrections to the 1 cm s⁻¹ level is in development, and we plan to present it in a future publication.

2.1.5. Convergence Criteria

Since the Doppler factor is the only strongly nonlinear parameter, the convergence criteria are established as the iteration requiring a δα_v × c smaller than 1 cm s⁻¹. Thanks to the simplicity of the model, only 3–6 iterations are typically required to reach convergence.

The typical nominal uncertainty of an echelle aperture with an S/N of 100 is between 2 and 15 m s⁻¹, strongly depending on the spectral type and aperture under consideration. Examples of template-matched observations for a G8.5V star (Tau Ceti) and an M4 dwarf (GJ 699) are shown in Figures 1 and 2.

**Figure 1.** Result of the fit to a small chunk of the spectrum of Tau Ceti (G8.5V). The average S/N of the observation is 180. This S/N would correspond to an rms of 0.55% if only photon noise were involved. The actual rms is 0.65%. Even though it is a very good fit, this slight excess of rms illustrates that assuming Poisson statistics is usually too optimistic to obtain the weight of each pixel and derive a realistic formal RV precision. A cubic flux normalization polynomial has been applied to match the blaze function of the observations to the blaze function of the template. The nominal RV uncertainty for this aperture (aperture 54) is 5.4 m s⁻¹.
Download figure:
Standard image High-resolution image

**Figure 2.** Results of the template-matching process for a small chunk of the spectrum of GJ 699 (M4V). The average S/N of the observation is 83, which should correspond to an rms of 1.22% compared with the actual 1.23% obtained. Even though the S/N is almost a third compared to the S/N of the G dwarf in Figure 1, the much more abundant Doppler information in the spectral features of the M dwarf gives a nominal uncertainty of 2.68 m s⁻¹ for this observation.
Download figure:
Standard image High-resolution image

2.2. Obtaining the Final RV Measurements

A final refinement is useful to produce the highest-quality RV. Because the template has uncertainties, there can be ambiguous RV zero-point offsets associated with each aperture. Also, some apertures might have extra uncertainty due to instrumental effects not accounted for by the formal uncertainties (extraction issues, problems with the wavelength solution, etc.). With the purpose of computing a more realistic weighted mean to each epoch, the zero point and nominal uncertainty of each order can be empirically reassessed as follows.

As a first step, the weighted mean μ_e of each epoch e is computed as

$\begin{eqnarray} \mu _e &=& \frac{1}{Z}\sum _{a=0}^{N_a}\frac{v_{e,a}}{\sigma _{e,a}^2}\,, \end{eqnarray} \tag{ 10 }$

$\begin{eqnarray} Z &=& \sum _{i=0}^{N_a} \frac{1}{\sigma _{e,a}^2}\,, \end{eqnarray} \tag{ 11 }$

where the subindex a runs over all the N_a apertures being analyzed, v_{e, a} is the RV measurement of the ath aperture at epoch e, and w_{e, a} is the corresponding formal uncertainty in the RV ( $\sigma _{\alpha _v}$ as obtained from Equation (9) multiplied by the speed of light). Z is the sum of all the weights and normalizes the weighted mean.

Let us now concentrate on one aperture, say, aperture 15. In each epoch e, we compute the difference between the epoch weighted mean μ_e and v_{e, 15}. The weighted average of this difference over all the epochs is the zero point of order 15, which we call v⁽⁰⁾₁₅. The weighted standard deviation of these differences over all the epochs is the new nominal uncertainty σ'_{e, 15}. This process is then repeated for all the orders to obtain the corresponding v⁽⁰⁾_i and σ'_{e, i}. After all the weights and zero points are obtained, the final RV measurement and its precision on a given epoch are given by

$\begin{eqnarray} \rm{RV}_e &=& \frac{1}{Z^\prime } \sum _{i=0}^{N_o}\frac{v_{e,i}-v^{(0)}_{i}}{\sigma ^{\prime 2}_{e,i}}\,, \end{eqnarray} \tag{ 12 }$

$\begin{eqnarray} \sigma _e^2 &=& \frac{1}{Z^\prime (N_o-1)} \sum _{i=0}^{N_o}\frac{\big(v_{e,i}-v^{(0)}_{i}\big)^2}{\sigma ^{\prime 2}_{e,i}}\,, \end{eqnarray} \tag{ 13 }$

$\begin{eqnarray} Z^\prime &=& \sum _{i=0}^{N_o} \frac{1}{\sigma ^{\prime 2}_{e,i}}\,. \end{eqnarray} \tag{ 14 }$

Note that, because v⁽⁰⁾_i and σ'_{e, i} are computed with respect to the initial guess of the epoch mean value, the zero-point correction does not remove offsets shared by all the apertures (i.e., Keplerian signals are not removed). In case the user is interested in the RV measurements without this re-normalization of the weights, the original μ_e and the individual RV measurements of each order are also provided as a data product by HARPS-TERRA.

As mentioned in the Introduction, these algorithms and all the code necessary to extract the HARPS spectra from the ESO reduced data products are included in the HARPS-TERRA software. HARPS-TERRA has been designed to achieve the highest precision while being computationally efficient. It takes between 1 and 5 s to process one full spectrum on a Linux operating system running on a 2.0 GHz processor. The construction of a high-S/N template (see the next section) takes around 5 minutes when using 200 spectra. A data set consisting of 100 spectra can be fully processed in less than 10 minutes.

2.3. Construction of the Template

In order to obtain maximal precision, one would need a template with the highest possible S/N. This can be achieved by carefully co-adding all the available spectra. Still, some precautions must be taken in the construction of this template. Below is an outline of the process we use to generate them.

In a first pass over all the spectra, the RV and flux normalization coefficients are obtained with respect to a preliminary template (highest-S/N observation). This preliminary RV measurement and the heliocentric motion of the observer are used to obtain the barycentric wavelengths of each observed spectrum at each epoch.

Spectra obtained in different epochs are not sampled at the same barycentric wavelengths and cannot be co-added without some kind of interpolation. To solve this, the preliminary template is used as a reference to generate a grid of regularly spaced reference wavelengths. The number of re-sampled reference wavelengths is four times the number of original pixels in each aperture of the preliminary template (this is 4096 × 4 = 16, 384). Each observed spectrum can then be interpolated on these reference wavelengths using a cubic spline. Finally, the flux normalization is applied and the template value at each re-sampled wavelength is computed as a 3σ clipped mean over all the epochs.

The flux normalization polynomial is only applied if the average S/N of a given aperture is higher than 5. S/Ns lower than 5 are not rare on the bluer apertures of M dwarfs, and the corresponding normalization polynomials are very unreliable. Because the blaze variability is of the order of 1%, it is safer to simply apply a scale factor and match the average flux of the observation to the average flux of the preliminary template.

Telluric features must also be removed from the co-adding. As in the RV measurement, the flux measurements coincident with telluric features deeper than 1% are masked out. Let us note that, if a star has been observed in different phases of the Earth motion around the Sun, this process allows us to generate a telluric free spectrum on regions with mild telluric contamination. When a significant number of spectra are available (N > 20), the result of this co-adding can be very spectacular, especially on M dwarfs whose spectrum has almost no continuum, and it is very hard to distinguish pure noise from real features (see Figure 2 as an example).

3. PERFORMANCE

Even though the described algorithms are relatively simple, the extraction of precise RVs still depends on a number of parameters that have to be tuned by hand. The two parameters we investigate here are (1) the optimal degree for the flux normalization polynomial and (2) the bluest echelle aperture to be used (e.g., M dwarfs can have 10–100 times more flux on the red than on the blue, and low-S/N spectra are typically more affected by systematic noise). HARPS-TERRA has been designed to be very flexible on all such parameters, but a general procedure to produce optimal results with minimal human intervention is still desirable. To illustrate the effect of these two parameters, we use observations of Tau Ceti (GJ 71, HD 10700, a very stable G8.5V dwarf), HD 85512 (GJ 370, a quiet K6V dwarf with a very low amplitude candidate planet), and Barnard star (GJ 699, a relatively quiet halo M4 dwarf). These three stars have abundant data in the ESO archive and bracket the highest-priority targets for the search of very low mass companions (G, K, and M dwarfs). The direct comparison of the final RVs obtained with CCF and HARPS-TERRA on a larger sample is given in Section 4. All the stars we discuss here are nearby and show a significant linear trend due to perspective acceleration (Zechmeister et al. 2009). In what follows, this perspective acceleration has been subtracted from the measured RVs. The basic parameters of all the stars discussed in the text are given in Table 2.

Table 2. Relevant Parameters of Each Star

Parameter	GJ 676A	Tau Ceti	HD 85512	Barnard's	Kapteyn's	Proxima	Eri	HD 69830
μ*_R.A.[mas yr⁻¹]^a	−260	−1721	461	−798	6505	−3776	−975	279
$\mu _{\rm Dec\rm{l.}}$ [mas yr⁻¹]^a	−184	854	−472	10328	−5731	766	19	−987
Parallax [mas]^a	61	274	90	548	256	772	311	80
V^b	9.58	3.5	7.651	9.51	8.85	11.05	3.73	5.95
K^b	5.82	1.79	...	4.52	5.05	4.38	1.78	4.16
Sp. type^b	M0V	G8.5V	K6V	M4V	M1V	M6V	K3V	G8V
Mass [M_☉]^c	0.71	0.78	0.69	0.16	0.27	0.12	0.82	0.86

Notes. A more detailed description can be found in the SIMBAD database and references therein. All quantities are given to the last significant digit. ^aProper motions and parallaxes from HIPPARCOS (van Leeuwen 2007) are required to subtract the perspective acceleration effect (see the text). μ*_R.A. corresponds to μ_R.A.cos δ and is the proper motion in the direction of increasing R.A. in a local tangent plane defined on the star's nominal coordinates at the catalog reference epoch. μ_R.A. is the obsolete coordinate-dependent definition of the secular change in R.A., which is singular at the celestial poles and does not have a direct physical interpretation. ^bUnless noted in the manuscript, V and K photometry and nominal spectral type are obtained from the SIMBAD database. ^cStellar masses for M dwarfs have been derived from Delfosse et al. (2000) using absolute K magnitudes. These masses should be accurate at the 5%–10% level. Masses of the K and G dwarfs are obtained from various references (see the section on each star).

Download table as: ASCII Typeset image

Tau Ceti. For this experiment we use 84 spectra of the 4000+ available in the HARPS-ESO archive. Since it is a very quiet and stable star, Tau Ceti had been used as a standard star by several HARPS programs. Not all such programs aim to achieve the highest RV precision. This results in a very heterogeneous sample of public observations with varying exposure times and a wide range of S/Ns. The S/N of some of the available spectra is so high (>500) that saturation or extraction problems seem to be seriously affecting the reddest orders. For a consistent analysis, we prepare a sub-sample from the HARPS high-precision survey for exoplanets in the southern hemisphere (e.g., Lovis et al. 2006; Mayor et al. 2009). This subset also belongs to the sample presented in Pepe et al. (2011) showing an rms of 0.92 m s⁻¹ over a time span of 5 years. Pepe et al. (2011) noted that, to achieve such precision, several spectra in a given night had to be averaged to mitigate the effects of stellar pulsation and granulation. Because this star is one of the 10 higher-priority targets in the HARPS-GTO Survey for Earth-like Planets, only the observations taken before 2009 are publicly available. In any case, we only need a statistically significant sample with consistently measured RVs using the CCF method as a reference. Still, to work with a more manageable sample, we only use the first observation of each night. The final sample contains exposure times between 30 and 150 s and a typical S/N in the range between 100 and 250 at 6000 Å. The CCF RVs for this sub-sample show an rms of 1.53 m s⁻¹. The rms from HARPS-TERRA using all the apertures is 1.52 m s⁻¹. Since we are only using one spectrum per night, this sub-sample shows, as expected, a larger rms than the nightly averages used in Pepe et al. (2011).

HD 85512 (GJ 370) is a quiet K6V dwarf with a very low amplitude candidate planet (Pepe et al. 2011) that could support liquid water if it were significantly covered by clouds (Kaltenegger et al. 2011). It is also a target of the current HARPS-GTO program to search for low-mass companions, so only the first 5 months of observations (2008 December to 2009 March) are publicly available. The typical integration times are between 400 and 600 s, and the S/N at 6000 Å ranges from 100 to 250. For our purpose here, it is only relevant that the CCF measurements show an rms as low as 1.10 m s⁻¹ using a K5 binary mask (a very stable star indeed). The HARPS-TERRA RVs from the same 122 observations show an rms of ∼1.0 m s⁻¹, which already represents a significant improvement. For comparison, the proposed planet candidate has an RV semi-amplitude of 0.8 m s⁻¹, and the detection is based on 250+ measurements in Pepe et al. (2011). A quadratic trend (also present in the Ca H+K S-index activity index) had to be removed to cleanly detect the candidate in Pepe et al. (2011). Neither trend nor planet signals were subtracted from the RVs discussed here.

Barnard's Star (GJ 699) is the star with the highest proper motion and the second-closest star system to the Sun. It is an M4 dwarf with halo kinematics (total velocity with respect to the Sun is ∼150 km s⁻¹), and it is known to be slightly metal poor. Even though it is a low-mass star and relatively faint in absolute terms, its proximity to the Sun (1.82 pc) allows us to obtain a typical S/N between 50 and 80 at 6000 Å in 900 s integrations. Twenty-two spectra are available in the HARPS-ESO archive over a time span of one year (2007 April to 2008 May). Even though it is classified as active (V2500 Oph), only occasional flares have been reported on it. Barnard's star was observed within the ESO-UVES Search of Low-mass Companions around M Dwarfs (Zechmeister et al. 2009). Those measurements demonstrated its RV stability down to 2.5 m s⁻¹ using the iodine cell technique. They also found that a possible RV wobble with a 40 day period was significantly correlated with the strength of its H_α emission. Therefore, we a priori expected that the RV measurements obtained with HARPS would also show some activity-induced jitter. The CCF RVs as extracted from the archive are obtained using an M2 binary mask and show an rms of 1.54 m s⁻¹. The RVs from the HARPS-TERRA RV measurements using the full spectrum have an rms of 1.23 m s⁻¹, which again represents a significant improvement. A quick-look analysis of the RVs does not show evidence of the previously reported periodicity at 40 days, but the number of observations is too low and the time sampling is too sparse to rule out signals in this period domain.

Two more stable M dwarfs. In addition to Barnard's star, we will briefly use HARPS RV measurements of two other M dwarfs: Proxima (GJ 551, M6Ve) and Kapteyn's star (GJ 191, M1V subdwarf). The data set on Proxima consists of 27 RVs taken between 2005 May and 2009 February. The CCF rms is 2.29 m s⁻¹, and recurrent flaring events can be detected in more than one epoch. The data set on GJ 191 contains 30 spectra taken between 2003 December and 2009 February. The CCF rms is 2.45 m s⁻¹. These two stars will be used to illustrate some features that seem to be common in M dwarfs other than Barnard's.

3.1. Flux Normalization Polynomial

Because RV measurements rely on the slopes of strong spectral features, the flux normalization correction is a key element for achieving the highest possible RV precision. For example, much better precision is obtained from deep sharp lines than from broad shallow ones. In a perfectly calibrated instrument, the blaze function should be constant over time. However, several instrumental/observational effects can cause variability of the effective blaze function. For example, due to atmospheric differential refraction, the photocenter of a star at the entrance fiber will depend on wavelength-causing wavelength-dependent flux losses as a function of the airmass. Also, bluer wavelengths are more efficiently dispersed by the atmosphere, adding additional airmass-dependent variability. In sum, the combination of several effects causes time-variable blaze function shapes at the level of a few percent, which (as shown below) can severely perturb the RV measurements.

Because we are interested in relative measurements, we only need to correct for the differential variations of the blaze as compared to the template. In Figure 3, we show the RV measurements of our three test stars when using different flux normalization polynomials. When only a constant normalization is applied (zero-degree polynomial), the scatter in the RVs of the G and K dwarfs is very large and has a clearly systematic behavior. The RV precision dramatically improves when using a linear flux correction (first-degree polynomial) on both stars. No significant improvement in precision is obtained by using polynomial degrees higher than 3. On the other hand, the RV measurements on the M dwarf are much less sensitive to the degree of the polynomial used. The improvement on the precision is modest but still significant with a linear correction, and optimal results are also obtained when a cubic polynomial is applied.

**Figure 3.** Radial velocity measurements as a function of time obtained on the three test stars using different flux normalization polynomials. Unless corrected, the flux variability strongly perturbs the RV measurements on the G star (left) as can be seen in the obtained rms when only an overall normalization factor is applied (order 0 in black). The effect of the flux variability is still significant for the K dwarf and has very little effect on the M dwarf. A first-order polynomial (blue) provides a significant improvement on both the G and K dwarfs, and a third-order polynomial (red) seems to reach a good compromise between computing efficiency and precision on all spectral types. RV offset has been added to the measurements of the zero- and first-degree polynomials to improve visualization.
Download figure:
Standard image High-resolution image

Sensitivity to the flux normalization correction as a function of spectral type was expected. The spectra of typical G and K dwarfs consist of well-isolated sharp lines against a smooth continuum. If not properly corrected, the changes in the slopes induced by the blaze variability on the continuum contain spurious Doppler information that strongly perturbs the Doppler measurement. On the other hand, the spectrum of an M dwarf is dominated by heavily blended molecular absorption bands that strongly dominate over the shifts induced by blaze variability. Given that high-degree flux corrections are computationally expensive, we set the nominal blaze correction to a cubic polynomial.

3.2. Most Useful Echelle Apertures

The final uncertainty in the RV measurement depends on many instrumental effects in addition to the formal statistical uncertainties and the stellar spectrum. K and, especially, M dwarfs have significantly less flux in the blue than in the red. If noise in bluer wavelengths was purely random, these apertures would be properly down-weighted through the formal uncertainties and the final measurement would be unaffected. However, low-S/N spectra can be more sensitive to instrumental systematic effects, so using arbitrarily low S/N data can be counterproductive. As an example, the M2 binary masks used by the HARPS DRS only use apertures redder than the 22nd one (λ > 4400 Å). Also, given that all the stars are active at some level and that activity should affect the measured RV differently at different wavelengths (Reiners et al. 2010), one could expect some apertures to provide more reliable measurements than others. As an example, star spots are known to have higher contrast at bluer wavelengths, so one would expect stronger RV stellar jitter on the blue. The wavelength dependence of the stellar jitter has been exploited in the past to rule out possible companions around young stars using complementary RV measurements in the near-infrared, e.g., see the RV measurements obtained in the H band by Huélamo et al. (2008) and Figueira et al. (2010) on TW Hya using the CRIRES spectrograph at the ESO/VLT. This section is devoted to developing a strategy to determine the wavelength dependence of the RV precision and defining the bluest aperture to be used for each spectral type.

Given that the typical formal uncertainty of a given aperture is never better than 2–3 m s⁻¹, it is hard to distinguish random noise from systematic effects (intrinsic to the star or instrumental) when looking at the RV measurements of single echelle apertures. Instead, we use the following procedure to assess the dependence of the precision with the bluest aperture used. First, we extract the RVs on all the apertures. Then we measure RVs on all the epochs as described in Section 2.2 but using only the apertures between 71 (the reddest available one) and some bluer one (say, aperture 68) and obtain the rms of these time series. We then repeat this process, adding one aperture at a time, and plot the obtained rms as a function of the bluest aperture used. This is illustrated in Figure 4 on our three test stars. On the right-hand side of each panel, the rms is expected to be higher because less apertures are used. However, if the noise were purely random, we would expect a reduced rms as more apertures are added to the blue. The K6V dwarf nicely follows this behavior.

**Figure 4.** Radial velocity rms as a function of the bluer aperture used. The rms derived from the public CCF measurements is indicated with a horizontal line. The bluer aperture used in the CCF analysis is marked as a thick vertical line. While the rms of the K6V dwarf decreases as bluer apertures are added, neither the M dwarf nor the G dwarf shows the same behavior. The origin of this rms excess in the blue is still unknown. The fact that it only shows on some very stable stars indicates that it is likely due to stellar-induced jitter rather than an instrumental effect.
Download figure:
Standard image High-resolution image

However, the M4V star behaves quite differently. The rms reaches a minimum (81 cm s⁻¹) at aperture 58 when only the 14 reddest apertures are used. Then, it starts increasing as more apertures are included toward the blue (∼1.2 m s⁻¹ at HARPS aperture 0). The typical S/N is still very significant on apertures bluer than 58 (S/N > 30), so instrumental systematic effects are not expected to be important at this level. This, combined with the fact that the K6V star does not show such an excess, suggests that a significant fraction of this extra noise is not instrumental but intrinsic to the star. For this star in particular, the 58th aperture (λ ∼ 6000 Å) seems to be the middle point where the activity-induced variability and the photon+instrumental noise equally contribute to the error budget. This behavior is also present in other stable M dwarfs. For example, we show the same diagrams for Proxima (M6V) and Kapteyn's star (M1.5V, GJ 191) in Figure 5. Because both stars are fainter than Barnard's, one has to use bluer apertures to reach the point where the instrumental noise meets this wavelength-dependent jitter. In an ideal situation the bluest aperture should be selected for each star, but if this effect is really related to stellar activity, this sweet spot is difficult to predict a priori. This behavior is under investigation and hints toward new ways of assessing if a signal could be activity-induced. Also, the wavelength dependence of the rms illustrates the importance of moving to redder wavelengths to efficiently search for low-mass companions around cool stars (Reiners et al. 2010). Figures 4 and 5 also show the rms obtained from the CCF measurements (horizontal line), illustrating that the template-matching technique works better even without "cherry picking" the bluer aperture to be used. As mentioned before, the CCF M2 binary masks used by the HARPS DRS do not consider apertures bluer than the 22nd. Given that aperture 22 also seems a reasonable choice for HARPS-TERRA, we will also use it as the default bluer aperture to be used when extracting RVs from M dwarfs.

**Figure 5.** rms as a function of the bluer aperture used on Proxima and Kapteyn's star (GJ 191) showing a behavior similar to Barnard's star. That is, the rms starts decreasing as bluer wavelengths are added until it starts increasing again (GJ 191) or does not change significantly (Proxima). For Proxima, the S/N below aperture 10 is so low (S/N <5) that adding such apertures to the analysis only adds systematic noise to the RV measurement. As a general rule, the cut at aperture 22 used by the HARPS DRS seems a reasonable assumption.
Download figure:
Standard image High-resolution image

Tau Ceti is known to be one of the most RV-stable G dwarfs. A bit unexpectedly, we find that the rms as a function of the bluer aperture also shows a minimum of 1.36 m s⁻¹ at aperture 40, and then it starts increasing again toward the blue (1.52 m s⁻¹ at aperture 0). Given that the S/N is still high (∼40 at order 0), again this should not happen in a perfectly stable star. As for the M dwarf, this indicates that some scatter is due to stellar activity and has a chromatic component (e.g., stellar pulsation should be achromatic). The rms obtained at the sweet spot (order 40th) can be used to estimate the magnitude of this chromatic noise as $1.36/\sqrt{2} \sim 0.96$ m s⁻¹, which is quite significant. For comparison, an Earth mass in the habitable zone on Tau Ceti would have an RV amplitude of 10–20 cm s⁻¹. Given that it is one of the most quiet Sun-like stars being surveyed for rocky planets, this issue requires a more detailed analysis, which is beyond the scope of this paper. Given the increased precision achievable by selecting the sweet-spot aperture, a reprocessing of the full data set in Pepe et al. (2011) with HARPS-TERRA could lead to a significant increase in the overall precision and enhanced sensitivity to lower-mass companions.

In conclusion to this section, we find that a cubic polynomial is sufficient to obtain an optimal flux normalization correction. The maximal precision on quiet stars (such as HD 85512) is obtained using the whole spectrum. However, in several stars, a moderate improvement in the overall precision could be achieved by finding the sweet spot between random noise (instrumental+photon) and the wavelength-dependent jitter. The determination of this sweet spot requires a better understanding of its physical origin and will not be further discussed here. To make a fair and consistent comparison, all the apertures will be used when analyzing spectra of G and K dwarfs. For M dwarfs, only apertures redder than aperture 22 (∼4400 Å) will be used. This bluest aperture used together with the cubic flux normalization polynomial discussed in Section 3.1 defines the standard setup of HARPS-TERRA when analyzing spectra from G, K, and M dwarf stars.

3.3. Ca ii H+K Activity Indicator

Pseudo-RV variations can be caused by stellar activity. As mentioned before, it is suspected that some apparent RV offsets are closely related to the magnetic activity of the star and related surface features such as spots. Also, enhanced magnetic activity can cause local or global changes in the convection patterns of the stellar surfaces. When convection is enhanced, hotter and bluer material emerging from the convection cells (e.g., stellar granulation) causes apparent blueshifts to the integrated stellar spectrum. For one reason or another, one could expect apparent RV jitter whenever the magnetic field of the star experiences changes. A detailed discussion on the topic can be found in Lovis et al. (2011) and references therein. Lovis et al. (2011) also provide the recipes for computing one of the most commonly used activity indicators, the Mount Wilson S-index, using HARPS data. This index measures the relative flux of the Ca ii H and K lines in emission (λ_H = 3933.664 Å and λ_K = 3968.470 Å) compared with a local continuum. These lines in emission are formed in the hot plasma of the chromospheres of stars, and, as a consequence, their intensity varies with the strength of the stellar magnetic field. We incorporated the automatic measurement of the S-index as one of the outputs of HARPS-TERRA. To obtain the correct flux estimates, it is necessary to know the absolute heliocentric RV of the star at a few hundred m s⁻¹ accuracy. Since HARPS-TERRA cannot provide this information, we use the median of the heliocentric CCF RV measurements to estimate the barycentric wavelengths of the lines on all the epochs. Because the H and K lines and the continuum bands defined in Lovis et al. (2011) appear in different echelle apertures (apertures 5 and 6), we use the re-sampled, blaze-corrected full spectrum also provided by the HARPS DRS. The full spectrum is then interpolated using a cubic spline, and the indices are computed on a regularly sampled grid of 0.01 Å conveniently Doppler-shifted to match the heliocentric RV of the star.

To validate this procedure, we obtained time series of the S-index on published HARPS stars, obtaining perfect agreement in all cases. The S-index measurements on Tau Ceti and HD 85512 are illustrated in Figure 6 and show the same behavior reported by Pepe et al. (2011). No coherent variability is observed on Tau Ceti (relative scatter is 0.7%), and long-term variability is observed on HD 85521. The S-index as a function of time is also plotted for Barnard's star. Barnard's (M4V) shows a similar amount of relative variability as the K dwarf except for a mild flaring event that doubles the Ca ii H+K flux in one of the epochs. In case a tentative signal is found, the time series of the S-index will be used to investigate possible correlation of the RVs with magnetic activity. An example of this is illustrated in Section 4.4, where we find that a promising long-period quadratic trend in the RV is correlated with an almost identical trend in the S-index. A popular form of the same activity indicator is the so-called R'_HK index. R'_HK represents the flux in emission of the Ca H+K lines over the total stellar spectrum and requires some assumptions on the spectral energy distribution of the star. Still, R'_HK is obtained as a linear relation from S, and, as a consequence, it will experience the same time variability. Since the computation of R'_HK requires additional information about the star (e.g., effective temperature and metallicity; Lovis et al. 2011) and does not provide extra information about the time variability, R'_HK is not provided by HARPS-TERRA.

**Figure 6.** S-Index time series on Tau Ceti, HD 88512, and Barnard's star. The mean value 〈S〉 and the relative variability of the index σ/〈S〉 are given for each star. One of the epochs for Barnard's star shows an S value two times higher than usual, which is indicative of a flaring event occurring during the observation. That point was excluded from the computation of the mean S value and its relative variability. The error bars represent the photon noise in the measurement of S.
Download figure:
Standard image High-resolution image

The HARPS DRS also provides three other activity indicators derived from the shape of the CCF. The first one is the so-called Bisector Span of the CCF (BIS). BIS is a measure of the asymmetry of the average spectral line and should correlate with the RV if the observed variability is caused by spots or plages rotating with the star (Queloz et al. 2001). The second index is the FWHM of the CCF and is a measure of the width on the mean spectral line. The variability of the FWHM is thought to be a direct consequence of changes in the convective patterns on the surface of a star, effectively changing the shapes of the integrated line profiles, but could also depend on other physical processes related to the stellar magnetic field. The third index is the contrast of the CCF (CONTRAST), which is sensitive to the changes in the depth of the average spectral line profiles. These three numbers provide important diagnostics to distinguish genuine Doppler signals from activity-induced periodicities. Since they are already computed by the HARPS DRS, they will not be further discussed here. Equivalent activity indicators optimally designed for the template-matching technique are in development and will be given in a future publication.

4. COMPARISON: CCF VERSUS TEMPLATE MATCHING

Comparing actual RVs is more informative than discussing which method (CCF or template matching) works better on theoretical grounds. This section is devoted to illustrating the performance of the template-matching approach compared with the pipelined HARPS-CCF-measured RVs. We discuss a few remarkable systems in terms of reported planet abundance but also their reported RV stability. We show that, while the performance with G and K stars is similar using both techniques, template matching works significantly better on M dwarfs and moderately active stars. Given the relative youth of the HARPS-TERRA code compared with the many years of refinement of the HARPS-CCF algorithms, we can only expect further improved precision in the future. In addition to alternative RV measurements, the template-matching approach should allow us to perform a new set of diagnostics such as the determination of the bluest aperture to be used for an optimal RV extraction already discussed in Section 3.2. The RV measurements relevant to each star are given in Tables 5 to 7.

Whenever signals are present and orbital fits are required, we use the SYSTEMIC interface (Meschiari & Laughlin 2010) as provided in 2010 August⁵ to obtain the orbital parameters and their uncertainties. SYSTEMIC allows the interactive adjustment of multi-planetary systems and is able to generate a large variety of data products and figures. When periodograms are required to illustrate the detection of a signal, we use a custom-made version of a least-squares periodogram. A sinusoidal signal is adjusted to a list of 10⁴ test periods between 1.1 and 10,000 days. The F-ratio statistic of the least-squares solution is then plotted against the period so the higher peak represents the most likely periodicity in the data. These periodograms are based on the definitions given in Cumming (2004), and they are formally equivalent to the so-called generalized Lomb–Scargle periodograms described by Zechmeister & Kürster (2009). The CCF and HARPS-TERRA RV measurements used in this section are given in the Appendix.

4.1. Proxima, a Flaring M4V

Proxima Centauri is the nearest stellar neighbor to the Sun and, therefore, has a special interest among the planet-hunting community and the public in general. So far, there has not been any serious claim of a companion around this star. It is the common proper-motion pair to the α Centauri binary, two stars with masses similar to the Sun in a long-period orbit. Proxima Cen has been intensively monitored by the ESO/UVES planet search program at a median precision of 3 m s⁻¹ using the iodine cell technique, and no significant periodicity was detected (Endl & Kürster 2008). Between 2004 and 2009, it has been observed 27 times with HARPS on different programs.

Visual inspection of Proxima's spectrum strongly suggests that it is an active star as several activity indicators in the HARPS wavelength range are in strong emission (H_α, NaD, Mg, Ca, etc.). Perhaps for this reason and because it is relatively faint, the star has not been monitored as intensively as other earlier-type nearby M dwarfs such as GJ 581 (Forveille et al. 2011a) or GJ 876 (Rivera et al. 2010). The HARPS-TERRA RVs have been obtained using the standard setup for M dwarfs. Figure 7 shows the CCF and HARPS-TERRA measurements as a function of time. On 2004 July 16, Proxima experienced an energetic event (probably a flare) causing all the activity indicators to go to strong emission. The flare event happened in the second epoch, and the emission lines (such as H_α) show intensities 2–10 times stronger than their quiescent state (S-index goes from a median value of 8.7 to 58). Still, neither the CCF nor the HARPS-TERRA measurement shows a significant offset on that particular epoch. Instead, the largest CCF offset (7.5 m s⁻¹) occurred 3 days later with no apparent counterpart in the activity indices. For comparison, the HARPS-TERRA measurement of this same epoch is only 3.2 m s⁻¹ away from the mean RV. When excluding the RV outlier, the rms of the CCF goes from 2.37 to 1.98 m s⁻¹ and TERRA goes from 2.05 to 1.88 m s⁻¹. The extra noise to be added in quadrature to the TERRA rms to match the rms of the CCF is 1.2 m s⁻¹ considering all the measurements, or 0.6 m s⁻¹ if the outlying event is excluded. An overview of the data sets on Proxima Cen and the other two M dwarfs discussed in the previous sections (Barnard's and Kapteyn's) is given in Table 3.

**Figure 7.** Measured radial velocities on Proxima using HARPS-TERRA (black) compared with those obtained using the CCF method (red).
Download figure:
Standard image High-resolution image

Table 3. Overview of the RV Measurements on the Three M Dwarfs without Planets Discussed in the Text

Parameter	Proxima	Barnard's	Kapteyn's
Mean S/N at 6100 Å	39	120	85
Number of RV measurements	27	30	22
First observation	2004 May	2003 Dec	2007 Apr
Last public observation	2009 Feb	2009 Feb	2008 May
rms_CCF [m s⁻¹]	2.37	2.45	1.51
rms_TERRA [m s⁻¹]	2.05	2.13	1.19
Extra noise in CCF [m s⁻¹]	+1.18	+1.21	+0.92
Optimal bluer aperture	37	29	56
rms at optimal aperture [m s⁻¹]	1.94	1.95	0.82

Download table as: ASCII Typeset image

Concerning possible signals, the periodogram of the HARPS-TERRA RVs shows a very marginal peak at 5.6 days that is also barely visible in the CCF values when the outlying RV measurement is removed. So, unfortunately, no promising signals are yet detected on Proxima Cen.

4.2. GJ 676A, an M1V with an Eccentric Gas Giant

In the previous sections, we have illustrated the improvement in the precision thanks to the use of the template-matching technique on stars with no planets. To illustrate that the increased precision is not an artifact of the data reduction process, we now analyze the spectra of a planet-hosting M dwarf with a large-amplitude signal. GJ 676A is orbited by a gas giant candidate with a period of 1060 days and an RV semi-amplitude of ∼120 m s⁻¹. GJ 676Ab was announced by Forveille et al. (2011b) and is one of the few gas giants detected around low-mass stars. The RV measurements provided by the HARPS DRS contain two epochs where the CCF failed to converge, giving spurious offsets of 38 km s⁻¹ and 76 km s⁻¹, respectively. The HARPS-TERRA measurements in these two epochs look perfectly reasonable. For comparison purposes, these two points were removed from the orbital analysis. As seen in Figure 8, the HARPS-TERRA RVs show essentially the same RV signal as the CCFs, demonstrating that our fitting procedure is not knocking out real RV offsets and that our rms estimates for stable stars are realistic.

**Figure 8.** Measured radial velocities using CCF (left) and HARPS-TERRA (right). Only the 34 epochs in the archive coincident with the published RVs are used for this fit. This figure clearly shows that the observed reduction in the rms of M dwarfs is not a spurious effect of the data analysis techniques we apply.
Download figure:
Standard image High-resolution image

After subtracting the best-fit solution for one planet, the CCF RVs show an rms of 6.31 m s⁻¹, while the HARPS-TERRA ones have an rms of 6.08 m s⁻¹. Given that the typical photon noise is of the order of 1–2 m s⁻¹, it is obvious that something else is happening with this star. As suggested by Forveille et al. (2011b), we added a linear trend to the fit and adjusted all the free parameters again. The improvement on both fits is quite significant (HARPS-TERRA rms is 3.13 m s⁻¹ and the rms from the CCF values is 3.77 m s⁻¹; see Table 4). As previously noted by Forveille et al. (2011b), the fact that the detected trend is significantly larger than the maximal expected acceleration due to GJ 676B strongly suggests the presence of an additional companion to GJ 676A with a period of at least several thousand days. Table 4 contains the best fit to the CCF and TERRA data sets. Only those epochs (34 of them) that are present in both data sets are included in both fits. The uncertainties correspond to the 68% confidence intervals as obtained using the Markov chain Monte Carlo (MCMC) algorithm included in SYSTEMIC. The MCMC jump lengths β of each parameter were tuned so the acceptance rate of all parameters was between 15% and 30% (Ford 2005), and 5×10⁶ MCMC iterations were used to generate the desired distributions and confidence levels in Table 4. The same data analysis technique will be used to characterize the uncertainties in the orbital parameters given in the forthcoming sections.

Table 4. GJ 676A Orbital Solution Using CCF and TERRA RVs

Parameter	CCF	TERRA
P [days]	1061.7 (2.3)	1060.2 (1.8)
K [m s⁻¹]	128.10 (66)	125.97 (58)
M₀ [deg]	210.13 (89)	208.77 (72)
e	0.326 (92)	0.331 (78)
ω[deg]	88.8 (1.5)	89.8 (1.2)
Linear trend [m s⁻¹ yr⁻¹]	8.48 (66)	8.99 (56)
Msin i [M_jup]	4.842 (25)	4.752 (22)
a [AU]	1.81	1.81

rms [m s⁻¹]	3.77	3.11
σ_{O − C}[m s⁻¹]^a	3.51	2.97
Extra noise^b	+2.13	...
N_obs	34	34
N_par	8	8
χ²	113.95	114.38
χ²/(N_obs − N_par)	4.38	4.39

Notes. The numbers in parentheses indicate the uncertainty in the last two significant digits of the parameter values. Uncertainties have been obtained using a Markov chain Monte Carlo approach (see the text). ^aWeighted rms of the residuals as computed by Pepe et al. (2011). ^bUncertainty that has to be added in quadrature to the CCF rms to match the rms of the TERRA measurements.

Download table as: ASCII Typeset image

Table 5. Differential RV and S-index Measurements for Tau Ceti, HD 85512, GJ 676A, and epsilon Eridani

JD	RV_TERRA	σ_TERRA	RV_CCF	σ_CCF	S-index	σ_S
(days)	(m s⁻¹)	(m s⁻¹)	(m s⁻¹)	(m s⁻¹)
Tau Ceti
2453280.550009	0.42615	0.51027	0.16668	0.30930	0.13932	0.00050
2453281.551756	1.58720	0.44181	1.69523	0.28281	0.14008	0.00045
2453282.551419	0.47161	0.52024	0.92778	0.33676	0.13800	0.00053
2453283.550154	1.18747	0.54740	0.46403	0.32926	0.13992	0.00052
2453284.793218	−0.16593	0.69264	−0.05192	0.60515	0.13798	0.00089

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Table 6. RV Measurements for Barnard's Star, Kapteyn's Star, and Proxima Cen

JD	RV_TERRA	σ_TERRA	RV_CCF	σ_CCF	RV_TERRA-Red	σ_TERRA-Red
(days)	(m s⁻¹)	(m s⁻¹)	(m s⁻¹)	(m s⁻¹)
Barnard's star
2454194.893890	0.00000	0.34651	0.57194	0.32376	0.00000	0.34651
2454196.883399	0.81911	0.35716	0.84845	0.30566	0.81911	0.35716
2454197.894815	−0.78946	0.33255	−0.38601	0.31239	−0.78946	0.33255
2454198.890257	1.50797	0.45543	2.70673	0.43272	1.50797	0.45543
2454199.918201	0.44577	0.32082	0.35208	0.26366	0.44577	0.32082

Notes. The last two columns are the RV measurements obtained using the optimal redmost part of the spectrum only.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Table 7. Measurements for HD 69830

JD	RV_TERRA	σ_TERRA	RV_TERRA-Red	σ_TERRA-Red	S-index	σ_S
(days)	(m s⁻¹)	(m s⁻¹)	(m s⁻¹)	(m s⁻¹)
2452939.874025	0.78514	0.44513	−0.44804	0.49607	0.14331	0.00049
2452943.876091	4.52766	0.54224	3.79835	0.50728	0.14377	0.00043
2452946.863070	1.66663	0.54780	−0.08756	0.62966	0.14381	0.00042
2452949.840269	5.06226	0.97553	1.66953	0.94835	0.13985	0.00070
2452949.842734	4.21729	0.92317	0.86545	0.95698	0.14056	0.00064

Notes. Instead of the CCF RVs, RVs using the redder part of the spectrum (see Section 4.4.1) are provided in this case.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

As we have seen on the other stable M stars, the post-fit rms from TERRA is significantly lower than the one obtained with CCF, giving further proof that, for M dwarfs, template matching significantly outperforms CCF. In this case, the uncertainty that has to be added in quadrature to the HARPS-TERRA measurements is 2.1 m s⁻¹, which is very significant. As noted in the discovery paper of GJ 676A, further low-mass planets could be present around GJ 676A in shorter periods. Given the significantly better accuracy, it is likely that such signals could already be detected by applying HARPS-TERRA to the complete data set. However, even though these results were already published, only 38 spectra are available through the database compared to the 69 used in Forveille et al. (2011b), so the full reanalysis cannot be done here. Also, we found that several RV measurements are missing in Forveille et al. (2011b). We checked that these missing measurements coincide most of the time with the outliers found in the archive. By comparison, no spectrum failed to provide a useful RV measurement using HARPS-TERRA, demonstrating the robustness of the algorithms.

4.3. Eridani, an Active K3V with a Gas Giant Companion?

epsilon Eridani is a K3V star of approximately 0.82 M_☉ (Butler et al. 2006) and a close neighbor to the Sun (3.2 pc). epsilon Eridani was first reported to be an RV variable by Campbell et al. (1988), and Walker et al. (1995) explicitly reported a possible variability with a period between 5 and 10 years. A 6.9 year periodicity was first reported in Cumming et al. (1999), but, at that time, it was considered suspicious given the high chromospheric activity on the star. Using additional RV measurements, Hatzes et al. (2000) proposed that the observed 6.9 year variability was most likely caused by the presence of a planet. The proposed planet would have a significantly eccentric orbit (e ∼ 0.7) and had a minimum mass of 0.86 M_jup. The reanalysis of several RV data sets combined with Fine Guidance Sensor/Hubble Space Telescope (FGS/HST) astrometry indicated that the companion was an actual planet with a true mass of ∼1.5 M_jup (Benedict et al. 2006). We want to remark that the rms of those early RV measurements was around 10–20 m s⁻¹ while the amplitude of the claimed candidate is about 18 m s⁻¹. Also, the astrometric amplitude reported in Benedict et al. (2006) is very close to the epoch-to-epoch systematic errors of FGS/HST. Because of the bias toward large amplitudes affecting astrometric measurements (Pourbaix 2001), this astrometric measurement should be understood as an upper limit. Evidence of gaps in the debris around epsilon Eridani suggests an additional very long period companion (not yet confirmed) that should have a negligible RV signature (e.g., Backman et al. 2009). epsilon Eridani is significantly more active than the Sun and has a strong stellar magnetic field (Rueedi et al. 1997). Magnetic-related activity (flaring, bright and dark spots, etc.) is the supposed source of the observed excess in RV variability in timescales of weeks. Stellar global magnetic cycles could also be responsible for part of the observed long-term variability.

The star was regularly observed by HARPS between 2004 and 2007. The archive contains 113 spectra of epsilon Eri. Ten spectra showed anomalous shapes, giving poor mismatches of the order of 20% in several echelle apertures. A closer look did show severe extraction problems on all of them. Three of them are from 2003 November 6, and three more are from 2004 February 4, all showing anomalous flux deficits between 5500 and 6000 Å. The other four are from 2004 September 13 and 2004 November 2, and all show very bumpy blaze function shapes probably related to extraction issues (according to the headers, these spectra were obtained on engineering time). After excluding these 10 measurements, the final sample we use contains 103 HARPS observations on 23 different nights.

The CCF measurements as extracted from the archive have an rms of 9.1 m s⁻¹, while the photon noise of each observation is of the order of ∼1 m s⁻¹. This variability shows no temporal coherence (no peaks in the periodogram), and therefore these observations seem to confirm that the stellar activity is contributing significantly to the observed scatter in the RVs. The spectra were processed using the standard HARPS-TERRA setup for K dwarfs. The resulting RV measurements show a remarkably lower rms (6.8 m s⁻¹) than the CCF RVs. In Figure 9, we present the nightly averages (23 equivalent epochs) overlapped to the nominal solution of Planet b as given by Benedict et al. (2006). In the case of the CCF, when the proposed planet is subtracted, the rms goes from 9.3 to 11.7 and the HARPS-TERRA rms increases from 6.8 to 10.0 m s⁻¹. Based on this, it appears that the planet is not real, or a major revision of its orbital solution needs to be obtained.

**Figure 9.** Twenty-three nightly averages of the 103 spectra used to measure the RVs. Red symbols represent the CCF measurements, and the HARPS-TERRA measurements are shown as solid black dots. No coherent RV variability coincident with previous orbital solutions is observed in either data set.
Download figure:
Standard image High-resolution image

Because of its proximity to the Sun and its brightness, epsilon Eridani has been observed by several programs with different instruments over the years. The six data sets available to date are provided in the current distribution of SYSTEMIC and were extracted from Benedict et al. (2006) and Butler et al. (2006). Let us note that, because the data used in Butler et al. (2006) were restricted to fewer measurements, the orbital solution presented there (only RV) was already quite different from the one reported by Hatzes et al. (2000) and Benedict et al. (2006). In order to check if there is still an orbital solution compatible with the new RV data, we add the new HARPS-TERRA measurements to a Keplerian fit with the other six data sets. In addition to the five Keplerian parameters (period, minimum mass, eccentricity, initial mean anomaly, and argument of the node), the model must include seven constants to account for the zero point of each instrument. We make a first tentative fit using the nominal reported uncertainties. The obtained solution is significantly different from the one reported by Hatzes et al. (2000), Benedict et al. (2006), or Butler et al. (2006). In particular, the new RV semi-amplitude is only ∼11 m s⁻¹ compared to the ∼18.0 m s⁻¹ given on all the previous studies and the orbital period changes from 6.9 to 7.25 yr. These changes are driven by the lack of coherent variability observed in the new HARPS-TERRA RVs. The obtained eccentricity is slightly smaller and the argument of the node is also forced to very different values (see Table 8). Especially in the most recent data sets with smaller formal uncertainties, it is obvious that the RV scatter is dominated by stellar noise rather than photon noise. Since the HARPS-TERRA measurements show the smaller rms, we use them to estimate the amount of stellar jitter that has to be added in quadrature to the nominal uncertainties in order to recover the obtained scatter. We find that this jitter amounts to ∼6.6 m s⁻¹. By adding this 6.6 m s⁻¹ in quadrature to all the reported uncertainties, we derive a more realistic orbital solution and corresponding uncertainties. The definitive orbital fit is quite similar to the one obtained without adding the jitter but has a χ² per degree of freedom of 1.15, indicating that 6.6 m s⁻¹ is a reasonable estimate for the stellar jitter. Figure 10 and Table 8 present the best fit to the seven RV data sets with the jitter included in the error bars. Even if this orbit seems a good fit to all the observations, we need to remark that the orbital solution is significantly different than the previously reported ones (e.g., the period P is at 4.2σ from previous estimates). We find this difference very suspicious and seem to indicate that the long-term RV variability of epsilon Eri is due to stellar activity cycles (non-strictly periodic) rather than a putative planet. Let us note that, even if the planet is real, the astrometric measurements of the orbital inclination and mass reported by Benedict et al. (2006) are no longer valid due to the much greater weight of the RV measurements in the determination of the orbit. In light of this result, and given that significant efforts are being devoted to attempt direct imaging of this candidate (e.g., Janson et al. 2007; Heinze et al. 2008), a reassessment of the allowed orbital parameters combining astrometry with RV measurements using modern Monte Carlo techniques (e.g., Reffert & Quirrenbach 2011; Anglada-Escudé et al. 2012b) is mostly needed. Also, a few additional HARPS observations fully covering the putative period should be sufficient to confirm that, at least, the observed variability is still present over a full orbital period. Precision RVs in the near-infrared over a full orbital period should provide definitive confirmation/refutation of the existence of a planet.

**Figure 10.** Best orbital solution satisfying all the radial velocities publicly available to date on Eridani. This solution is significantly different from the one presented in previous works. The new HARPS-TERRA measurements are the black dots on the right. The naming of the data sets follows the prescriptions given in Benedict et al. (2006).
Download figure:
Standard image High-resolution image

Table 8. Orbital Solution for epsilon Eridani Including Seven Data Sets

Parameter	CCF
P [days]	2651 (36)
K [m s⁻¹]	11.8 (1.1)
M₀ [deg]	09 (12)
e^a	0.40 (11)
ω[deg]	141.4 (9.8)
Msin i [M_jup]	0.645 (58)
a [AU]	3.51

rms [m s⁻¹]	12.6
σ_{O − C}[m s⁻¹]^b	10.7
N_obs	359
N_par	12
χ²	401.02
χ²/(N_obs − N_par)	1.15

Notes. A stellar jitter of 6.6 m s⁻¹ was added in quadrature to all nominal uncertainties. Parameter values represent the least-squares solution to the last two significant digits. The numbers in parentheses represent the uncertainty in the last two significant digits. The uncertainties represent the 68% confidence level intervals as obtained from a Bayesian MCMC. ^aAsymmetric distribution. 0.20 < e < 0.68 with a 99% confidence level. ^bWeighted rms of the residuals as computed by Pepe et al. (2011).

Download table as: ASCII Typeset image

In summary, the precision in the RV obtained using the TERRA code provides a significant reduction on the measured rms compared to the public CCF RVs. Given that epsilon Eridani is an active star, the stellar lines used in the CCF might be peculiar when compared with quiet K dwarfs, contributing to the excess of scatter. Even though stellar activity is the most likely cause for most of the observed jitter, we do not see the chromatic jitter effect seen on M dwarfs or Tau Ceti. Figure 11 shows how the rms always decreases as bluer apertures are included. This indicates that the leading effect causing the observed RV jitter on active stars is quite different from the one affecting quieter dwarfs.

**Figure 11.** rms vs. bluest aperture used in the RV measurement. Even though it is an active star, the rms always decreases as bluer apertures are included, indicating that the chromatic component of the rms is not a dominant source of RV jitter. The CCF measurements also use all the apertures.
Download figure:
Standard image High-resolution image

4.4. HD 69830, Three Neptunes and a Long-period Trend?

HD 69830 made the news in 2006 (Lovis et al. 2006), being the first planetary system hosting three Neptune-mass planets. HD 69830 is a G8V star slightly less luminous and less massive than the Sun (0.6 L_☉ and 0.86 M_☉, respectively). The detection was based on 74 HARPS epochs taken between 2003 October and 2006 January. The final fit to the solution had an rms below 1.0 m s⁻¹.

The ESO archive contains 529 spectra of the star. The headers of two spectra were corrupted, rendering those two observations unusable (both observations are from 2007 February 5). The resulting 527 spectra were taken between 2003 October and 2008 March. The S/N at 6000 Å ranges from 100 to 300, and the exposure time varies between 180 and 400 s. The spectra were processed using the standard HARPS-TERRA setup for G dwarfs. By comparison, the CCF RV shows several outliers with RV offsets of several km s⁻¹ away from the average (nine spectra on three different nights). When removed, the rms of the full CCF sample is still rather high compared with that of the TERRA sample (∼6.74 m s⁻¹). By direct inspection of the CCF RV measurements, we could see that several RVs show negative offsets of the order of 20 m s⁻¹. We found that all such measurements (67 of them) correspond to spectra processed by the HARPS DRS with a G2 binary mask instead of the K5 binary mask used for the majority of the observations (454 of them). Because the HARPS DRS is not public, we could not consistently reobtain all the CCF RVs. We show all the RV measurements obtained with HARPS-TERRA and the CCF ones obtained using the K5 binary mask in Figure 12.

**Figure 12.** Five hundred twenty-seven radial velocity measurements obtained with HARPS-TERRA (black) compared with the 454 RVs obtained using a K5 binary mask with the CCF method (red circles). The CCF RVs are shifted 15 m s⁻¹ to improve visualization.
Download figure:
Standard image High-resolution image

The purpose of this section is to illustrate that HARPS-TERRA can detect signals at the limit of the demonstrated HARPS stability. Given the general good agreement among the two data sets and for the sake of simplicity, only the HARPS-TERRA measurements will be discussed below.

In order to perform the orbital analysis, we consolidate these 527 observations into 176 nightly averages. The periodogram of these RVs directly shows the three signals reported by Lovis et al. (2006). The best-fit solution to these three planets compares well to the solution given in the discovery paper and has an rms of 1.1 m s⁻¹.

However, the periodogram to the residuals to the three-planet fit still shows a strong peak at 490 days (see Figure 14) and an even higher power beyond 1000 days. By visual inspection of the residuals (see Figure 13), one can see that the long-period power is due to an apparent parabolic shape in the RV residuals. When we tried solving for the 490 day signal, we found that the coverage in the phase of this solution was very spotty. This is characteristic of long-period signals aliased with the seasonal availability of the star (see Dawson & Fabrycky 2010 for a detailed discussion of the yearly alias). Given that the power at longer periods is higher than the 490 day peak, from now on we assume that the long-period trend is the most likely signal left in the data.

**Figure 13.** Left: RV residuals to a fully Keplerian fit to the planets already reported by Lovis et al. (2006). A significant quadratic trend (red line) is observed in the data. As reported by Lovis et al. (2006), the early measurements (left) have larger uncertainties due to lower S/N and other instrumental issues. Right: S-index measurements as obtained from the original 529 spectra. A clear quadratic trend is also observed on this index pointing toward activity as the origin of the observed long-term variability in the RV.
Download figure:
Standard image High-resolution image

To assess the odds of this signal being generated by random noise, we compute its empirical false alarm probability (or FAP) as follows. We generate 10⁴ synthetic data sets by keeping the observed epochs but doing random permutations of the residuals to the three-planet fit. While these synthetic data sets will have the same distribution of random errors, the random permutations destroy the temporal coherence of any signal present. We then compute the periodogram for each synthetic data set and count the number of times we obtain a power higher than the one we find in the real data. This only happened once, indicating that the FAP is low (∼0.01 %), which means that this trend is statistically significant and cannot be ignored. This method of computing FAP is described in more detail in Cumming (2004).

In order to assess the possible cause of this trend, we examined the measurements on the S-index as described in Section 3.3. Unfortunately, the S-index shows a trend with the same shape and relative variability as the RV curve. This is characteristic of spurious offsets caused by the magnetic activity cycle of the star and has been observed in other stable stars such as Tau Ceti (Pepe et al. 2011). A tentative fit of the RV residuals provides a period of ∼8500 ± 2000, which roughly matches the expected duration of the activity cycle of a quiet star similar to our Sun. Even though the signal could be the combination of a long-period planet and the activity-induced offset, we cannot conclude that there is solid evidence for an additional companion to HD 69830. Once a circular orbit is subtracted from the data, the signal at 490 days completely disappears from the periodogram, confirming that both peaks in Figure 14 correspond to the same aliased trend.

**Figure 14.** Periodogram to the residuals of HD 69830 to a three-planet fit (top) and to a four-planet fit (bottom). Even though there is a suggestive periodicity around 490 days, this signal is most likely an alias of the long-period parabolic trend seen in Figure 13. When the best circular orbit to the long-period trend is removed, no additional signals are present in the data.
Download figure:
Standard image High-resolution image

The number of observations used here is significantly larger than that available at the time of discovery of the first three candidates. Also, only the long time span of the observations (∼5 years) allows us to detect the quadratic trend and its correlation with the S-index variability. As shown in the bottom panel of Figure 14, no more periodicities can be inferred from the residuals for the four-signal solution. Removing the trend with a circular orbit leaves an rms 0.92 m s⁻¹. Note that, even with the trend, this is one of the most RV-stable G dwarfs observed by HARPS. As a final note, we analyzed the CCF RV obtained with the K5 binary mask using the same procedure. We recovered the same orbital solution for the three reported candidates, but both the long-period trend and its alias at 490 days appear with less significance. While one would be tempted to claim that HARPS-TERRA also obtains higher precision in this case, this extreme cannot be confirmed here for three reasons: (1) the CCF data set contains fewer measurements, (2) the rms of both measurements after removing the first three planets is almost identical, and (3) such precision is at the limit of the long-term stability of HARPS and any signal at this level (especially a long-period signal) has to be taken with due caution. An updated orbital solution is given in Table 9, and the RV measurements phase folded to the period of each planet candidate are shown in Figure 15.

**Figure 15.** Best fit to the four signals present in the HD 69830 measurements obtained with HARPS-TERRA. The first three planets were already reported by Lovis et al. (2006). A circular orbit for Planet e has been assumed to generate this plot. The rms of the residuals to the four-signal fit is only 0.92 m s⁻¹.
Download figure:
Standard image High-resolution image

Table 9. Orbital Solution for the Three Planet Candidates plus Trend Detected on the RVs of HD 69830 When Using All the Echelle Apertures

Parameter	b	c	d	(Activity?)
P [days]	8.6687 (12)	31.645 (28)	202.2 (1.6)	8500^a
K [m s⁻¹]	3.46 (11)	2.61 (22)	1.85 (13)	12.35^a
M₀ [deg]	228 (10)	118 (15)	39 (17)	144^a
e	0.06^b (05)	0.08^b (05)	0.190^c (76)	0 (fixed)
ω[deg]	24 (23)	192 (40)	172 (32)	0 (fixed)
Msin i [M_jup]	0.0315 (11)	0.0366 (22)	0.0473 (35)	1.12^a
Msin i [M_⊕]	10.00(32)	11.60 (69)	15.0 (1.1)	355.4^a
a [AU]	0.079	0.186	0.641	7.7 AU^a

rms [m s⁻¹]	0.92
σ_{O − C}[m s⁻¹]^d	0.88
N_obs	176
N_par	19
χ²	120.35
$\frac{\chi ^2}{N_{\rm obs}-N_{\rm par}}$	0.76

Notes. The proposed parameter values are obtained from the χ² solution. The numbers in parentheses correspond to the uncertainty in the last two significant digits in the parameter values (68% confidence level intervals). Statistical quantities at the bottom correspond to a circular orbital fit to the quadratic trend. All the orbital elements are referred to in the first epoch of observation at JD₀ = 2452939.87402 days. ^aOrbital values of the equivalent circular orbit. Note that the most likely explanation for this signal is the RV offset induced by the magnetic cycle of the star. ^bCompatible with circular orbit. ^cSlightly significant eccentricity. ^dWeighted rms of the residuals as computed by Pepe et al. (2011).

Download table as: ASCII Typeset image

4.4.1. Wavelength Dependence of the Signals in HD 69830

We performed an additional test to assess the reality of the aforementioned trend. Given the similar spectral type of HD 69830 (G8V) to Tau Ceti (G8.5V), we asked ourselves if the chromatic jitter effect detected in the rms of M dwarfs and Tau Ceti could be exploited to investigate the nature of the trend. To do this, we first plotted the rms as a function of the bluer aperture used and found its minimum (see Figure 16). We note that, although there are three low-amplitude candidates contributing to the initial rms, Keplerian signals should be achromatic and, therefore, a minimum in the rms should still be present when the wavelength-dependent noise is added in quadrature. Figure 16 shows that such a minimum in the rms is found at aperture 37, which is similar to the optimal bluer aperture found for Tau Ceti.

**Figure 16.** Radial velocity rms as a function of the bluer aperture used on HD 69830. Even though this raw rms contains the signal of three low-amplitude companions, the contribution of the Keplerian signals and the chromatic jitter should add in quadrature, providing an optimal bluer aperture to be used.
Download figure:
Standard image High-resolution image

**Figure 17.** Periodogram to the residuals of HD 69830 to a three-planet fit (top) when only apertures redder than 37 are used to produce the RV measurements. No hint of a trend or any other periodicity is seen in the data.
Download figure:
Standard image High-resolution image

We then used the RVs obtained considering only apertures redder than 37 and derived the full Keplerian solution for the first three planet candidates. The obtained orbits were compatible with the previously reported orbits. However, when we computed the periodogram of the residuals to the three-planet fit, we found that the quadratic trend and the corresponding 490 day alias were completely gone. Also, the rms to the three-planet fit was already very low (1.01 m s⁻¹), clearly indicating that the secular signal was mainly driven by apparent RV offsets at the bluest wavelengths. Even though instrumental effects might be involved, we think this result strongly suggests that (1) this chromatic jitter is intrinsic to the star, (2) it is a significant source of the stellar RV noise observed in other relatively quiet stars, and (3) it might pose a fundamental limit to the maximum RV precision achievable on G stars.

5. SUMMARY AND CONCLUSIONS

We present a new method for obtaining precision RVs from wavelength-calibrated public HARPS spectra. Our new velocities compare well with those obtained with the CCF method and are able to detect the same signals reported by previous HARPS discoveries. For stable G and K dwarfs, the difference between CCF and HARPS-TERRA is at the level of ∼0.5 m s⁻¹ rms. Still, the template-matching approach seems to be less sensitive to offsets induced by the stellar activity (e.g., see Section 4.3 on epsilon Eridani) and requires almost no assumptions on the nature of the star. Because the template-matching technique makes more optimal use of the Doppler information in the stellar spectrum, HARPS-TERRA provides a significant increase in accuracy when applied to the heavily blended spectra of M dwarfs.

We have shown the importance of correcting for blaze function variability (flux normalization polynomial) to achieve sub m s⁻¹ precision. While the CCF method needs some preprocessing of the spectra to account for such variability, least-squares template matching performs such a correction in a self-consistent way.

We find that several stars show excess variability when the bluer echelle apertures are used. Even if the S/N is typically lower in the blue, this increase in jitter should not happen with perfectly stable stars. We find significant evidence that this chromatic jitter is likely related to the star itself rather than to an instrumental effect and that the wavelength dependence of the RV offsets can be exploited to confirm or rule out suspicious candidate signals. Further research is necessary to assess the nature of this excess.

We also have demonstrated that HARPS-TERRA can detect already reported signals, and we have applied it to a number of interesting stars with abundant public data. While we can reproduce other detections well, the new HARPS-TERRA measurements do not firmly confirm the planet candidate reported around epsilon Eridani. The orbital solution allowed by the new observations is significantly different compared with the previously reported ones, casting some doubts on the reality of this candidate.

We also report the detection of a quadratic trend in the residuals of HD 69830. Despite the high significance of the signal, we find that it correlates well with a similar trend in the S-index. We also find that this trend completely disappears when only the redder half of the spectrum is used to derive the RV measurements. This leads us to conclude that the most likely explanation for the observed trend is the activity cycle of HD 69830 inducing wavelength-dependent RV shifts that are stronger toward the blue. Although we do not report a new planet, this example demonstrates that HARPS-TERRA is also able to robustly detect and diagnose unreported signals at the limit of the HARPS instrumental stability. Given that HD 69830 is a nearby star, a planet candidate with such a long period would lie around 0 farcs 7 from the central star and might be imaged with the next generation of adaptive optics systems (e.g., Lagrange et al. 2010). Even though everything points to stellar activity as the most likely explanation for the trend, we provide an estimate of the equivalent circular orbit just in case.

We have shown that the template matching on stabilized spectrometers requires few assumptions compared to the elaborate binary masks required by the CCF method. The capability of reproducing precision RV measurements with two different data analysis methods is also a powerful ally for double-checking the significance of very low amplitude signals. Given the significant increase in precision achieved on low-mass stars, it is likely that HARPS-TERRA can uncover undetected low-amplitude signals in already-existing data sets (Anglada-Escudé et al. 2012a).

The research of G.A. has been supported by the Carnegie Fellowship Postdoctoral Program. We thank Mathias Zechmeister, Steve Vogt, Alan Boss, and Alycia Weinberger for useful discussions. R.P.B. gratefully acknowledges support from NASA OSS Grant NNX07AR40G and from the Carnegie Institution of Washington.

APPENDIX: RADIAL VELOCITY MEASUREMENTS

This Appendix contains tables with the most relevant time series used in the paper. Unless stated otherwise, the first two columns are the RVs as obtained with HARPS-TERRA using the standard setup for each spectral type. The corresponding CCF values are also provided for comparison purposes, if necessary. An offset (average RV) has been subtracted from all the RVs to improve readability, and all RV measurements are Doppler offsets measured in the solar system barycentric system. The perspective acceleration effect has also been removed from all the presented RVs. No nightly averages have been applied. The tables also contain the S-index values in the Mount Wilson system as measured by HARPS-TERRA. In a few cases, the RV measurements using the redder part of the spectrum are provided instead of the S-index or the CCF measurements. Check each table caption for further information. RV measurements used in Anglada-Escudé et al. (2012a) to report the planetary system around GJ 667C are also included in Table 10.

Table 10. Differential RV Measurements on GJ 667C Used in Anglada-Escudé et al. (2012a)

JD	RV_TERRA	σ_TERRA
(days)	(m s⁻¹)	(m s⁻¹)
2453158.764365	−3.130	1.490
2453201.586793	−11.436	1.596
2453511.798845	−7.158	1.506
2453520.781048	−3.902	1.661
2453783.863347	0.630	1.210

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

THE HARPS-TERRA PROJECT. I. DESCRIPTION OF THE ALGORITHMS, PERFORMANCE, AND NEW MEASUREMENTS ON A FEW REMARKABLE STARS OBSERVED BY HARPS

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION