© ESO 2018

1. Introduction

The radial-velocity (RV) technique was the first method that allowed the detection of exoplanets orbiting solar-type stars. By measuring the stellar RV variations induced by the gravitational pull of orbiting planets, this method allows to indirectly measure the mass of those companions. After the first detection of 51 Pegasi b reported in Mayor & Queloz (1995), the field of exoplanets boomed and new planetary detections, with smaller and smaller masses were announced (e.g. Butler et al. 1999; Santos et al. 2004; Lovis et al. 2006; Queloz et al. 2009; Pepe et al. 2011; Anglada-Escudé et al. 2016; Feng et al. 2017).

With the improvement of high-resolution spectrographs over time, and with more performant techniques to extract the RV information from raw frames, the RV precision reached a limit of 1 m s⁻¹, or perhaps slightly better for spectrographs like HARPS (Mayor et al. 2003) and HARPS-N (Cosentino et al. 2012). At this precision, we start to characterise spurious RV stellar signals induced by physical phenomena happening on the surface of solar-like stars.

On the timescale of a few minutes, stellar oscillations modify the velocity of the stellar surface and therefore induce a spurious RV signal (e.g. Kjeldsen & Bedding 1995; Bouchy & Carrier 2002; Arentoft et al. 2008; Bazot et al. 2012). On timescales from minutes to days, we can see the effect of granulation, which is driven by the variation of convection on the stellar surface (e.g. Del Moro et al. 2004; Del Moro 2004; Lefebvre et al. 2008; Dumusque et al. 2011a; Cegla et al. 2013; Meunier et al. 2017). On a timescale similar to the stellar rotation period, stellar magnetic activity, responsible for the appearance of inhomogeneities like spots and faculae on the stellar surface induces a semi-periodic RV signal that can mimic a planetary signal (e.g. Saar & Donahue 1997; Queloz et al. 2001; Desort et al. 2007; Meunier et al. 2010; Dumusque et al. 2011b, 2014; Borgniet et al. 2015). Finally, on the timescale of several years, solar-like magnetic cycle induces long-term drift in RV measurements (Lindegren & Dravins 2003; Dumusque et al. 2011c; Lovis et al. 2011a).

Those stellar signals perturb the detection of tiny planetary signals. We recall here that the Earth induces a RV variation of 0.1 m s⁻¹ on our Sun, which is an order of magnitude smaller that the known sources of stellar signals described in the preceding paragraph. Over the years, different techniques have been developed to try mitigating the impact of those stellar signals. For example, it is common to observe stars for 15 min to allow averaging out the signal induced by stellar oscillations (e.g. Pepe et al. 2011). Taking several measurements per night of the same target helps mitigating the impact of granulation (e.g. Dumusque et al. 2011a). Probing the variation observed in activity indicators like the $\log \left(R_{\mathrm{HK}}^{\prime }\right)$ (Noyes et al. 1984), the BIS SPAN (Queloz et al. 2001) or H α (e.g. Bonfils et al. 2007; Robertson et al. 2014) helps in disentangling stellar activity from planetary signals. In addition, using statistical techniques like Gaussian Processes (e.g. Haywood et al. 2014; Rajpaul et al. 2015; Jones et al. 2017; Delisle et al. 2018) or Moving Average (e.g. Tuomi et al. 2013) can mitigate to a certain level the impact of the red noise induced by stellar activity. Finally fitting the RVs with the long-term trend observed in the $\log \left(R_{\mathrm{HK}}^{\prime }\right)$ activity index reduces the impact of solar-like magnetic cycle on the RVs (e.g. Dumusque et al. 2011c, 2012; Delisle et al. 2018). Although these different techniques have been applied on several RV data set and have been able to detect interesting planetary systems, it seems that nowadays we are facing the obstacle of stellar activity, which is extremely difficult to overcome. We therefore need to rethink the way the RVs are derived if we want one day to mitigate stellar activity at the 0.1 m s⁻¹ precision level that should be reached by the ESPRESSO spectrograph (Pepe et al. 2014).

In this paper, we investigate a new way of deriving RVs. The idea behind this study is to use the wealth of information contained in high-resolution stellar spectra rather than averaging all the information into the few parameters derived from the cross-correlation functions (Baranne et al. 1996). Some preliminary works from Davis et al. (2017), Thompson et al. (2017) and Wise et al. (2018) show that stellar activity affects spectral lines in a different way. This is not surprising as each spectral line have its how sensitivity to temperature, to magnetic field strength and to convection velocity, three physical parameters that are strongly modified in the presence of stellar activity. We therefore investigate in this paper the possibility of measuring the RV on each individual spectral line.

In Sect. 2, we discuss the different techniques that have been used so far to derive precise RVs and discuss their advantages and limitations. In Sect. 3 we describe the technique we used in this paper to measure the RV on each individual line, and in Sect. 4 we compare our new RVs with the ones derived using the HARPS Data Reduction Sofware (DRS), which is the gold standard for deriving precise RVs from HARPS spectra. In Sect. 5, we use our new RV extraction procedure to investigate night-to-night RV offsets induced by the non-stability of HARPS wavelength solutions. In Sect. 6, we demonstrate that stellar activity affects each spectral line in a different way, and that by using a smart selection of spectral lines when measuring the RVs, it is possible to either boost the RV stellar activity signal by a factor of 2, or mitigate it by a factor of 1.6. Finally, we discuss our results and conclude in Sect. 7.

2. Measuring precise stellar radial velocities

In the last 15 years, we could reach the radial-velocity precision of 1 m s⁻¹. Two different types of high-resolution spectrographs were developed and improved to reach this goal:

The first type of instrument consists of slit spectrographs not well controlled in temperature and pressure, which implies that atmospheric changing conditions induce hundreds of m s⁻¹ variations due to the modification of the instrument point spread function (PSF). To be able to reach the m s⁻¹ precision, the stellar light passes through an iodine cell, that imprint the absorption spectrum of iodine on top of the stellar spectrum. The iodine emission spectrum is affected by PSF variations in the same way as the stellar spectrum and therefore serves as a reference scale to measure precise RVs. By using a complex reduction, fitting simultaneously (i) a stellar spectrum deconvolved from the instrument PSF, (ii) an extremely high-resolution spectrum of the iodine cell, normally obtained with an Fourier Transform Spectrograph, and (iii) the PSF of the spectrograph, it is possible to measure on small chunks of the stellar spectrum, generally a few Angström, the local RV shift. Finally averaging the RV shifts of all the different chunks of the stellar spectrum gives stellar RVs with precision close to the m s⁻¹ (e.g. Butler et al. 1996).

The second type of instrument consists of fibre-fed spectrographs put in vacuum (Baranne et al. 1996), which allows for an extreme stability in temperature and pressure. The PSF of the instrument is stabilised, and therefore very small RV variations are observed, generally bellow the m s⁻¹ level. To correct for these small instrumental variations, the spectrum of a calibration lamp, thorium-argon (Th-Ar) or Fabry-Perot (FP) étalon, is simultaneously recorded on the detector with the stellar spectrum. To reach a precise RV measurement, several techniques have been used like cross correlation (Baranne et al. 1996; Pepe et al. 2002b), maximum likelihood template matching (e.g. Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015) or least square deconvolution (Donati et al. 1997). Those methods are used to average out the RV signal of all the spectral lines together, which allows to reach the m s⁻¹ RV precision.

Deriving precise radial velocities using stabilised RV instruments

The first reference of using the numerical cross-correlation technique to measure precise RVs in the context of exoplanets can be found in the paper describing the ELODIE spectrograph (Baranne et al. 1996). Once the stellar spectrum is recorded on the CCD, a cross-correlation is performed using a synthetic template (Baranne et al. 1996). The result of this transformation, namely the cross-correlation function (CCF), is an average line profile that looks like an inverse Gaussian. The RV of the star is derived by fitting an inverse Gaussian to the CCF and taking its mean. We note that this cross-correlation technique is inherited from the method used before CCDs exited. During those times, a physical template with holes at the position of every stellar spectral line was shifted in the focal plane of the spectrograph to construct physically the cross-correlation function that was then recorded using a photometer (i.e. for the CORAVEL instrument, Baranne et al. 1979)

The numerical cross-correlation technique is still used nowadays on stabilised spectrograph like HARPS, HARPS-N and soon ESPRESSO. The technique was improved in Pepe et al. (2002a), were the authors show that rather than using a binary mask, 0 for the continuum and 1 for the centre of spectral lines, weighting the mask by the depth of the spectral lines improve the RV precision. HARPS and HARPS-N proved that the technique was able to deliver in a simple and robust way sub- m s⁻¹ RV precision on stellar spectra. In addition to simplicity, the cross-correlation with a numerical mask has the advantage of directly giving a precise RV estimate for a given star without requiring any preceding observation.

It was however demonstrated that for M dwarfs, the template matching technique gives better results than the CCF technique using a weighted synthetic mask optimised for M dwarfs (e.g. Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015). In this technique, the RV is derived by performing a maximum likelihood estimation between each stellar spectrum and a high signal-to-noise ratio (S/N) master spectrum built from all the available stellar observations, in which the Doppler shift is a free parameter. For M dwarfs, the template matching technique gives more precise RVs because due to the low temperature of these stars, spectral lines are everywhere in the spectrum and therefore a master stellar spectrum can access to all the possible RV information present in the spectrum, what a synthetic mask with a limited number of lines cannot do. This approach gives however a worse RV precision when applied to G and K dwarfs. This is probably because in this case, the template matching technique is sensitive to the noise present in the continuum where no RV information is present. Note however that this could be solved by masking spectrum continuum regions in the master stellar spectrum.

An advantage of the cross-correlation technique is that we can measure, in addition to the RV, the full width at half maximum (FWHM) and the bisector of the CCF. Several studies have shown that those indicators are useful to assess the stellar or planetary nature of a RV signal (e.g. Queloz et al. 2001; Desort et al. 2007; Meunier et al. 2010; Dumusque et al. 2014). Indeed, a planet only induces a pure Doppler-shift of the stellar spectra, without changing the shape of spectral lines, and therefore without changing the FWHM and bisector of the CCF, which is not the case for all known stellar signals.

By using a template that contains all the available lines in a stellar spectrum, the CCF or template matching techniques can reach an extreme RV precision. However, by averaging all the spectral lines together, any strong effect that specific spectral lines might have will be diluted in the final product, and thus impossible to observe. The problem is therefore not linked to the CCF or template matching techniques, but by the fact that the RV information of all lines are averaged out together. For example, in Dumusque et al. (2015) the authors discovered that several stars observed with HARPS were presenting a RV signal of a few m s⁻¹ with a period of one year, in phase with the Earth barycentric RV. By looking closer at the RV of each spectral line, they were able to show that only a few spectral lines, the ones crossing the HARPS detector stitchings, were affected by a yearly signal with an amplitude up to a hundred m s⁻¹. By averaging out the signal of those lines with all the others not affected, the residual signals was of the order of a few m s⁻¹.

3. Measuring the radial velocity of individual spectral lines

The new approach presented in this paper to derive precise RVs is an intermediate solution between measuring locally the RV like in the iodine cell technique, to avoid diluting the RV signal of every spectral line, and using a stellar master spectrum without considering the stellar continuum, which might be the limitation of the template matching technique on G–K dwarfs. The idea is to measure the RV on each spectral line using the template matching technique, and then to perform a weighted average over all those lines to obtain a precise RV measurement for a given star.

The process of measuring the RV on each spectral line, starting from HARPS reduced 2d stellar spectra, can be split in different steps that we will describe further:

• for a given star, correct each HARPS 2d stellar spectrum from known systematics,

• build a high S/N 2d stellar master spectrum by co-adding these corrected 2d stellar spectra,

• select the spectral lines for which we want to measure their RVs,

• measure the RVs on those individual spectral lines,

• and finally combine the RV information of all spectral lines to get a precise RV measurement.

3.1. Spectra corrections from known systematics

To get a stellar continuum that is flat, we must correct the stellar spectrum from the blaze of the instrument. On HARPS, the daily calibrations measure the blaze using a tungsten lamp and the information is recorded in a HARPS.XXX_blaze_FIBRE.fits, XXX being the time of the calibration and FIBRE being A or B depending if the file is referring to the blaze of the science or reference fibres, respectively file. The name of this file is recorded in the header of every stellar observation, under the keyword ESO DRS BLAZE FILE. To correct for the effect of the blaze, we simply have to divide the stellar spectrum from it.

After dividing by the blaze, the stellar continuum is much flatter, however, we still see small linear trends on each spectral order. This is because the stellar energy distribution is not corrected for. Taking into account the effective temperature of the star we are studying, we measure, using the black-body radiation, the flux that is recorded by each pixel on the CCD and correct for this effect. We note that to perform this correction, we need to know the wavelength of each pixel on the CCD, which is given by the wavelength solution. This solution is measured by the HARPS afternoon calibrations and is recorded in the header of each stellar observation. A total of 288 polynomial coefficients are saved under the keywords ESO DRS CAL TH COEFF LLX, where X goes from 0 to 287. The wavelength λ_i,j of pixel j in order i can be obtained using the following formula:

$${\lambda }_{i,j}=P_{4*i}+P_{4*i+1}\times j+P_{4*i+2}\times j^2+P_{4*i+3}\times j^3,$$(1)

where P_X corresponds to the value of the keyword ESO DRS CAL TH COEFF LLX, i varying from 0 to 71, and j from 0 to 4095.

The relative flux between the blue and red parts of the spectrum of a given star, commonly referred to as a colour index, is correlated with atmospheric extinction that varies with airmass and changing atmospheric conditions (e.g. Bourrier & Hébrard 2014). The airmass effect is corrected for before taking the measurement using an atmospheric dispersion corrector, however the effect due to changing weather conditions cannot be accounted for a priori. This change in relative flux will put different weights as a function of time on the RV of each spectral order and thus can induce a correlation between the RV measured and colour index. To avoid this colour dependence, the total flux in each order is divided by the flux of a reference spectrum of similar spectral type, and the residuals are fitted with a 5th order polynomial, which coefficients are saved in the header of the observed spectrum under the keywords DRS FLUX CORR COEFFX, X going from 0 to 4. To get a similar colour index for each spectrum of a given star, and therefore obtain RVs that are insensitive to colour, we divide each spectrum with the corresponding 5th order polynomial.

3.2. Build a high S/N master stellar spectrum

Once each HARPS 2d spectrum is corrected from known systematics (see Sect. 3.1), building a high S/N stellar master 2d spectrum is a rather easy task, however some small details need to be taken into account. HARPS 2d spectra are recorded in the reference frame of the observatory, which velocity with respect to a given star changes as a function of time due to Earth orbiting the Sun and rotating on its own axis. This implies that the stellar spectra of a given star moves on the HARPS CCD as a function of time. To get the spectral lines of all spectra overlapping each other, we need to change the reference frame of the observation to the solar system barycentre. This can be done by Doppler shifting every stellar spectrum using the barycentric Earth RV (BERV) measured by the HARPS data reduction software (DRS) and saved in the header of each observation under the keyword ESO DRS BERV. Spectra can also be Doppler shifted by the presence of companions orbiting a given star. Therefore, in addition to correct for the BERV, we also correct for the RV measured by the HARPS DRS. Note however that since the width of a pixel on the HARPS CCD is equivalent to 820 m s⁻¹ only massive close-in binaries and hot Jupiters will have a significant impact on the high S/N stellar master spectrum.

By red or blue-shifting the spectra to correct for the BERV and the effect of potential companions, the minimum and maximum wavelength of each order will be different for each observation. We therefore search for the maximum of the minimum wavelength and the minimum of the maximum wavelength for each order over all observations, and set those values as being the minimum and maximum wavelength of each order for the high S/N stellar master spectrum.

Each stellar spectrum will have its own wavelength solution, however, only one wavelength solution can be used for the stellar master spectrum. We therefore choose the wavelength solution of the spectrum having the minimum BERV, and then linearly interpolated all the other spectra to this fixed wavelength solution. Choosing the wavelength solution for which the BERV of the star is minimum will assure that the median RV of all the spectral lines will be the closest to 0.

Even after correcting the spectra from known systematics (see Sect. 3.1), some orders still show small drifts. To remove those drifts before co-adding all the spectra together to build the high S/N stellar master, we fit the continuum of each order with a linear drift and correct for it. For reasons that we did not investigate, but probably due to a bad cosmic ray correction, the slope over an order can reach unusual values in some cases. To prevent adding bad spectra or systematics in the master stellar spectrum, we reject all spectra for which at least one slope over the orders 30–71 is further away than 5-σ from the mean of all slopes. We do not reject spectra with bad slopes for orders smaller than 30 as fitting for the stellar continuum in the blue part of the spectrum is challenging.

To have a high S/N stellar master spectrum that does not present significant instruments systematics, we also rejected spectra for which the S/N measured in order 10 (∼4040 Å) is smaller than ten and for which the S/N measured in order 60 (∼6110 Å) is larger than 450. This prevents of considering spectra that could be affect by the non-linearity of the detector but also by inefficiency in the charge transfert when the S/N is very low. On HARPS, the atmospheric dispersion corrector corrects for airmass colour effect up to an airmass of 2. To be conservative and to prevent adding to the stellar master template spectra that could be contaminated by atmospheric effects, we select only observations with an airmass smaller than 1.5. Finally, to get a master template that is not affected by stellar activity, we select only spectra for which the calcium activity index is smaller than $\log \left(R_{\mathrm{HK}}^{\prime }\right)<-4.95$. Those cuts in airmass and $\log \left(R_{\mathrm{HK}}^{\prime }\right)$ are very restrictive and we choose them here as we will analyse, in Sect. 4, stars that have thousands of spectra available. We will see in Sect. 3.4 that the noise in the master spectrum will be neglected compared to the noise in each individual spectrum. Therefore, as a rough estimate, at least 25 spectra must be selected to build the master as this will imply that on average the noise of the master will be five times smaller than the noise in each individual spectra. For some stars, the restrictive cutoff in airmass and activity level selected here will never be reached. In those cases, we would advise to select at least the 25 less active spectra with airmass smaller than 2. This is however something that should be studied further.

After doing all those corrections and rejecting spectra that could exhibits some systematics, we can build the high S/N stellar master spectrum by co-adding all the individual corrected spectra. The obtained stellar master spectrum is oversampled by a factor 82 using a cubic spline to reach a resolution of 10 m s⁻¹. We did not investigate what is the optimal oversampling factor to be used. But for computation performance, we preferred doing a dense interpolation at this stage using a cubic spline and then a linear interpolation for the final refinement to obtain a master spectrum which is on the same wavelength scale as each stellar spectrum (see Sect. 3.4). Because this stellar master spectrum is at very high S/N, we will be able, in Sect. 3.4, to neglect the noise added during this oversampling process.

3.3. Select the spectral lines to derive RVs

Three different synthetic templates are used by the HARPS DRS to measure the RV with the CCF technique: a G2, a K5 and a M2 mask. Those templates have been optimised and do not consider spectral lines strongly affected by magnetic activity, like the Ca II H and K, Hα or Mg II lines, but also spectral lines in the red part of the spectrum strongly contaminated by tellurics. When measuring the RVs of G or K dwarfs, we will always derive the RVs on all the individual spectral lines present in the K5 synthetic template to account for the maximum number of possible lines. If no line is present at the position defined by the mask due to a higher stellar effective temperature, the RV precision measured on this line will be extremely bad and the line will be rejected at a later stage (see Sect. 3.5).

3.4. Measure the RV on individual spectral lines

Let’s consider that we want to measure, RV_i, the RV of spectral line i with central wavelength λ_i. We first need to define a window around λ_i on which we will measure the RV of this line. To prevent windows overlapping in the case of close spectral lines, mainly in the blue part of the spectrum, but also to select the maximum number of pixel to reach the best precision in RV, we selected a window width of 16 pixels, which corresponds to 0.17 and 0.30 Å in the bluest and reddest parts of HARPS spectra, respectively. For simplicity, we selected only a fixed number of pixels for all spectral lines. Note however that the window width could be optimised further as spectral lines in the red are much less packed than in the blue.

The second step consists in selecting the high S/N stellar master spectrum chunk S_{ref, i}(λ) centred on λ_i and that present a width of 16 pixels. We record then the minimum ν_min and maximum wavelength ν_mac of S_{ref, i}(λ).

Then for each spectrum j with its own wavelength solution, we select the spectrum chunk S_{j, i}(λ) between ν_min and ν_mac and perform linear interpolation on the edges so that S_{j, i}(λ) for all times j present the same boundaries as S_{ref, i}(λ).

To calculate the RV shift of S_{j, i}(λ) with respect to the reference S_{ref, i}(λ), we need an estimate of the noise present in S_{j, i}(λ). This noise has two components that should be added in quadrature. One that is due to photon noise and that is equal to the square root of the flux, and another one that comes from the CCD read-out and dark current, estimated to 12 photo-electrons on HARPS, which can become significant at low S/N. To calculate the photon noise in HARPS spectra, we need to be careful to use the raw flux before correcting for the blaze and the stellar spectral energy distribution (see Sect. 3.1). Otherwise, we will put too much weight on spectral lines that are on the edges of the orders and that have lower flux due to the blaze. To have the correct flux and photon noise, we adjust the flux of S_{j, i}(λ) so that it is equivalent to the raw uncorrected spectrum.

To compute RV_{i, j}, the RV of line i at time j, we use the template matching method described in Bouchy et al. (2001), which uses the change in flux between two spectra and the derivative of one of them to measure the RV shift. This method can be applied if the two spectra are on the same wavelength scale, if the shift between the two spectra is smaller than the sampling, and if the reference spectrum on which we calculate the derivative has a high S/N so that we can neglect the noise in this spectrum and in its derivative. In our case, we apply this method on S_{j, i}(λ) and S_{ref, i}(λ), the latter being used as the reference. Following Eq. (2) of Sect. 2 in Bouchy et al. (2001), we can show that:

$$S_{j,\mathrm{ }i}(\lambda )\mathrm{ }=A\left[S_{\mathrm{ref},\mathrm{ }i}(\lambda )+\frac{\partial S_{\mathrm{ref},\mathrm{ }i}(\lambda )}{\partial \lambda }\delta \lambda \right]$$(2)

where A accounts for the difference in flux between S_{j, i}(λ) and S_{ref, i}(λ). We note that here we use spectra as a function of wavelength and not pixel. Moreover, we could use a more complex model to fit a zero-point offset to account for possible background contaminations, and/or a slope. We tried with an offset and this was not significantly improving our results, so we preferred to remove this additional parameter. This is however something that should be investigated when analysing the RVs of fainter targets that the ones we analyse here. Regarding the slope, we do not expect any significant change from spectrum to spectrum on such short wavelength windows.

Due to the BERV and potential orbiting companions, the shift between S_{j, i}(λ) and S_{ref, i}(λ) is often more than the sampling, which is 820 m s⁻¹ in the case of HARPS. To use the technique described above, we need to reduce this shift to a level smaller than the sampling, which is done by applying a Doppler shift on S_{ref, i}(λ) to add to the reference the RV shift induced by the BERV and potential orbital companions.

The template matching method requires that the template S_{ref, i}(λ) is on the same wavelength scale as the spectrum chunk S_{j, i}(λ) of one observation. This is done by estimating the template at the wavelength of the observation using linear interpolation¹. Once this is done, we fit Eq. (2) with a linear least-square algorithm taking into account the noise in S_{j, i}(λ), to get the best solution for A and A δλ. The RV shift between S_{j, i}(λ) and S_{ref, i}(λ), RV_{i, j}, is then derived using the Doppler shift formula:

$$\frac{\delta V}{c}=\frac{\delta \lambda }{\lambda }\Rightarrow {\mathrm{RV}}_{i,\mathrm{ }j}=\frac{c}{\lambda }\frac{A\delta \lambda }{A},$$(3)

where c is the speed of light. The square root of the diagonal elements of the covariance matrix obtained by the linear least square gives us the errors for A and A δλ, σ_A and σ_{A δλ}, respectively. The error for RV_{i, j} is obtained from propagating these errors:

$${\sigma }_{{\mathrm{RV}}_{i,j}}=\frac{c}{\lambda }\sqrt{{\left[\frac{1}{A}\right]}^2{\sigma }_{A\delta \lambda }^2+{\left[-\frac{A\delta \lambda }{A^2}\right]}^2{\sigma }_A^2}\cdot$$(4)

For all the stellar spectra, we estimate RV_{i, j} and σ_{RVi, j}, and therefore get the RV of spectral line i with its error as a function of time. Because the high S/N template spectrum is used as a reference, and because this spectrum is an average of all the stellar spectra corrected from the BERV and from the RV contribution of orbiting companions (see Sect. 3.2), the RV of each individual line will be centred around zero.

3.5. Combining the RV information of all spectral lines to get the stellar RV

In the preceding section, we explained how we could get the RV of each individual line. Those RVs can be useful to study stellar physics at the level of each spectral line. The RV precision on each line is however not very good, and reaches in the best case ∼10 m s⁻¹ (see Fig. 5) when analysing the spectra of the extremely bright target α Cen B (V = 1.33). To reach the m s⁻¹ precision, we need to combine the RV of all the spectral lines together.

To combine the RV information of all spectral lines, we first must perform some data cleansing. Due to imperfections on the CCD, saturation from cosmic rays that are not always corrected properly, spectral lines from our K5 line list (see Sect. 3.3) that do not appear in a given star due to a different spectral type and probably other unknown effects, the RV measured on each spectral line can show some significant outliers. To reject bad RV measurements and bad spectral lines, we first remove from the RV of each line the effect of the BERV and the effect of potential companions by subtracting the RV measured by the HARPS DRS. We perform then the following cleansing:

• We perform a 6-σ clipping on the χ² of all the spectral lines in one observation obtained when fitting Eq. (2) using a linear least square (see Sect. 3.4). We repeat the process twice, which eliminated spectral lines badly fitted.

• We perform a 6-σ clipping on the RVs of all the spectral lines in one observation. We however keep all the RVs for which the value plus six times its error bar is compatible with the median RV of all the lines, which is close to 0 by construction. This removes RV measurements that are too far from the median of all lines in one spectrum. We repeat the process twice.

• We fit on the RVs of each spectral a second order polynomial as a function of time, and perform a 6-σ clipping, without rejecting points for which the error bar is compatible with the polynomial fit at 6-σ.

• We reject all spectral lines for which more than 1% of the data are rejected by the three preceding cleansing.

• We reject the 0.3% of spectral lines with the worse RV precision.

• We remove all lines for which the ratio between the standard deviation of the RVs and the median RV error, what we call further the quality factor, is larger than two. Those lines have a significant signal, which is not expected as the RV precision on a line is poor and we expect stellar signals to be at a lower level.

• We finally reject spectral lines for which their median RV is not compatible at 3-σ with the distribution of the mode of all lines. We use the mode here as the right and left sides of the distribution can have different tails.

After curating the RVs of all spectral lines using the criteria defined above, we combine the RV information of all the remaining lines to get a precise measurement of the stellar RV at time j. This is done by performing a weighted average of the RVs measured on each spectral line, putting as weight their respective RV errors:

$${\text{RV}}_j=\frac{{\displaystyle {\sum }_i\left[\frac{1}{{\sigma }_{{\text{RV}}_{i,j}}^2}{\text{RV}}_{i,j}\right]}}{{\displaystyle {\sum }_i\left[\frac{1}{{\sigma }_{{\text{RV}}_{i,j}}^2}\right]}}·$$(5)

4. Comparison of the RVs derived from individual lines with the HARPS DRS RVs

In this section, we compare for three stars the RVs derived from our new method, that is analysing the RVs of each spectral line and then combining their information to get the RV of the star, with the RVs derived from the HARPS DRS using the CCF technique.

4.1. Radial velocities of τ Ceti

The star τ Ceti (HD 10700) has been observed since the beginning of HARPS in 2003. Over 15 years, more than 10 000 spectra have been obtained, with a very dense sampling. This G8 dwarf has the particularity of being inactive, with an activity level that always stayed below $\log \left(R_{\mathrm{HK}}^{\prime }\right)=-4.94$ over 15 years. This dwarf is probably in a similar state that the Sun was during its Maunder Minimum, which corresponds to the period around 1645–1715 during which sunspots on its surface became exceedingly rare. τ Ceti is therefore the perfect target to compare our new RVs with the RVs derived from the HARPS DRS.

Before deriving the RVs with our new RV extraction procedure, we have to build the master stellar spectrum that will be used as a reference to measure the RV of each single spectral line as a function of time (see Sect. 3.2). For τ Ceti, we selected spectra for which the S/N in order 10 was larger than ten, for which the airmass was below 1.2 and for which the activity level was below $\log \left(R_{\mathrm{HK}}^{\prime }\right)=-4.95$. We therefore considered 5743 spectra to build the master spectrum, representing 50.8% of all available spectra. We note that out of all those spectra, 15% of them are coming from the five-day asteroseismology run performed in October 2004 (Teixeira et al. 2009). Although not done here, we should investigate if selecting a lot of spectra taken at the same time and therefore at the same particular state of the star does not have an impact on the final RVs. After rejecting the spectra for which the slope of the continuum in some spectral orders was too high (see Sect. 3.2), we were left with 3769 spectra (33.4%), that were used to build the master stellar spectrum of τ Ceti.

After computing the RV for each spectral line using the formalism described in Sect. 3.4, we cleaned the data from any outlier before combining the RV of all the spectral lines as explained in Sect. 3.5. In total, out of the 7331 lines available in the K5 line list used, 701 lines were rejected in this cleansing process, thus slightly less than 10%². Rejecting only 10% of the lines might seems rather small, however, using the HARPS K5 template line list as a starting point in our analysis allows to remove most of the lines that would be rejected at this level (see Sect. 3.3).

Finally, before analysing the newly derived RVs, we have to correct for the drift of the instrument, as explained in Appendix A. As we can see in this section, our new RV extraction procedure on calibration spectra gives extremely similar results than the HARPS DRS, therefore we used our new drift measurements to correct for instrumental systematics in the τ Ceti RV data set.

In the top plot of Fig. 1, we compare the RVs for τ Ceti measured using our new RV extraction procedure with the RVs estimated by the HARPS DRS. The middle plot illustrate the difference between these two sets of RVs. The change of the HARPS fibres from circular to octagonal, which happened on JD = 2457174 (1st June 2015), did not have the same impact on the two reductions. However, besides this offset, the RVs are extremely similar with a rms for the RV difference of 0.40 m s⁻¹ when considering the data before the fibre change. We could remove this offset by building two master spectra, using either spectra taken before or after the fibre change, by measuring the RV of all lines in each spectrum using the corresponding master spectrum and finally by measuring the offset between the two masters. Looking, in the bottom plot of Fig. 1, at the periodograms of the two different sets of RVs for the data obtained before the fibre change, we see once again a similar behaviour. It seems however that there is slightly more power in our new RVs at periods between a hundred and a thousand days. Although we speak of signals with amplitudes of the order of 0.5 m s⁻¹, we still investigate the cause of this difference in Appendix B. In resume, it seems that our new RV extraction procedure is slightly more sensitive to the effect induced by lines crossing CCD stitchings (Dumusque et al. 2015) and by micro-tellurics (e.g. Cunha et al. 2014).

Fig. 1. Top panel: comparison of the RVs of τ Ceti derived using our new RV extraction procedure (red) with the RVs estimated from the HARPS DRS (green). The red vertical dashed line corresponds to the change from circular to octagonal fibres that was performed on HARPS on the 1st of June 2015 (JD = 2457174). Middle panel: difference between the RVs derived using our new RV extraction procedure and the RVs from the HARPS DRS. We show in the legend the standard deviation of the RV difference for the measurements preceding the fibre change. Bottom panel: periodograms of the RVs derived using our new RV extraction procedure (red) and the RVs derived by the HARPS DRS (green). To calculate those periodograms, we rejected the data after the fibre change to prevent low frequency signals induced by the observed offset. The two periodograms have been normalised in such a way that their FAP of 1% coincide.

Even though we see small differences in the periodogram of the RVs derived with our new RV extraction procedure and the HARPS DRS, when looking at the histogram of the RVs in Fig. 2, we see very similar Gaussian-like distributions with a similar mean and standard deviation, confirming once again that the two set of RVs are extremely similar.

Fig. 2. Histogram of the RVs of τ Ceti, HD 128621 and HD 10180 derived using our new RV extraction procedure (red) and the RVs derived using the HARPS DRS (green). The legend shows the standard deviation obtained for both data sets.

4.2. Radial velocities of α Cen B

Another star that has been intensively observed with HARPS at extreme precision is α Cen B (HD 128621). Contrary to τ Ceti, the HARPS data of α Cen B show a solar like magnetic cycle, with periods of minimum and maximum of activity. Those data are therefore important to check how our reduction compares with the HARPS DRS in the case of RVs affected by stellar activity.

In Dumusque et al. (2012), the authors show that the data of α Cen B are sometimes, depending on the observing conditions, contaminated by the light of α Cen A which is a few arcseconds away on the sky. We used here the same criterion as in Dumusque et al. (2012) to reject contaminated measurements. In addition, it was also shown in this paper that the data of α Cen B exhibit a significant signal at one year and its different harmonics due to CCD stitching (Dumusque et al. 2015). To remove those perturbing signals, we rejected from our new RV extraction procedure spectral lines that were closer to CCD stitchings than 48 pixels. We did the same for the RVs derived by the HARPS DRS, by using a cross-correlation template that had those lines removed (Dumusque et al. 2015). We note that in Bauer et al. (2015) and Coffinet et al. (2018), the authors develop a new recipe to derive HARPS wavelength solutions based on the physical size of the HARPS CCD pixels, which solves this problem of stitchings.

To build the master of α Cen B, we selected spectra for which the S/N in order 10 was larger than 100, for which the airmass was below 1.5 and for which the activity level was below $\log \left(R_{\mathrm{HK}}^{\prime }\right)=-4.95$. In addition, we rejected also spectra contaminated by α Cen A (see Dumusque et al. 2012), which gave us in the end 1200 spectra, which corresponds to 5, 7% of all available spectra, that were used to build the master stellar spectrum.

During the data cleaning process described in Sect. 3.5, a total of 310 spectral lines were rejected out of 7317. In addition, 1071 lines were not considered when deriving the final RVs because they were too close to CCD stitchings.

We compare in Fig. 3, as for τ Ceti, the RVs obtained with our new RV extraction procedure and with the HARPS DRS. Once again, we see extremely similar RVs, with a RV difference standard deviation of 0.43 m s⁻¹. The periodogram of both reduction are extremely similar, with this time signals at one year and its first harmonic that are slightly less significant when using our new RV extraction procedure. The histogram of the RVs, shown in Fig. 2, show once again very similar distributions. Therefore, in the case of stellar activity, our new RV extraction procedure gives extremely similar results than the HARPS DRS.

Fig. 3. Same as Fig. 1 but for the RVs of α Cen B.

4.3. Radial velocities of HD 10180

After testing our new RV extraction procedure on α Cen B, a star that shows solar-like activity, and on τ Ceti, a very inactive star, we test here if our reduction can be used to look for planetary signals. Although 5 planets have been announced orbiting τ Ceti (Feng et al. 2017), most of the signals found are really at the detection limit, and thus two different reductions of the data might give different results. We therefore analyse here a simpler case, the HD 10180 planetary system which is composed of 6 Neptune-like planets (Lovis et al. 2011b).

To build the stellar master spectrum used as reference to measure the RVs of each spectral line, we selected spectra that have a S/N in order 10 larger than 10, for which the airmass is smaller than 1.5 and for which the activity level is smaller than $\log \left(R_{\mathrm{HK}}^{\prime }\right)=-4.95$. After rejecting bad spectra as explained in Sect. 3.2, we are left with 210 spectra that are used to build the master stellar spectrum, which represents 62.5% of all spectra available.

During the data cleaning process described in Sect. 3.5, a total of 497 spectral lines were rejected out of 7317.

In Fig. 4, we compare the RVs derived using our new RV extraction procedure with the RVs from the HARPS DRS. Like in the case of α Cen B and τ Ceti, we find very similar RVs, with a standard deviation of the RV difference equals to 0.57 m s⁻¹. The periodograms of the two set of RVs looks very similar, likewise the histogram of the RVs shown in Fig. 2.

Fig. 4. Same as Fig. 1 but for the RVs of HD 10180.

To check if our new RV extraction procedure gives as precise orbital parameters as the HARPS DRS for the planetary signals present in the data, we first rejected measurements obtained after the fibre change, as we estimated that nine data points is not enough to robustly fit for the offset induced by this intervention. We also did not consider, when fitting for the different signals, the Earth-mass planet announced in Lovis et al. (2011b) at 1.2 days, as this discovery is at the limit of detection. After selecting the same data and model for both the RVs derived with our new RV extraction procedure and with the HARPS DRS, we used the MCMC tool available on the DACE platform³ to derive precise orbital parameters for the planets orbiting HD 10180. We used the orbital parameters derived in Lovis et al. (2011b) as initial guess for our MCMC fit.

The results of our orbital parameter estimation are summarised in Table 1. As we can see, all the orbital parameters derived using both RV data sets are compatible within 1-σ. The log Posterior is also compatible within 1-σ, and the RV residuals after removing the best fitting model only differs by 0.02 m s⁻¹. Therefore, our new RV extraction procedure gives RVs that are as good as the HARPS DRS RVs to look for planetary signals.

5. Night-to-night RV offsets due to the wavelength solution

As explained in Sect. 3.1, the HARPS DRS derives the wavelength solution by fitting on each spectral order a third order polynomial on the known wavelength of the Th-Ar emission lines (Redman et al. 2014). On the edges of each order, the flux on the Th-Ar emission spectrum is low due to the blaze, and thus the polynomial is not well constrained at those locations and can vary significantly from night-to-night. Therefore, some stellar spectral lines that appear on the edges of the spectrograph orders can have a wavelength solution that varies of several dozens of m s⁻¹ from night-to-night. Still due to the blaze, those spectral lines will have a small flux, and will therefore not contribute strongly to the HARPS DRS RVs, which combines the RV information of all spectral lines together. Nevertheless, the HARPS DRS RVs should still be affected by night-to-night RV offsets below the m s⁻¹ level.

In Appendix C we use our new RV extraction procedure to demonstrate that the RVs of spectral lines falling on the edges of the HARPS detector are indeed affected by dozens of m s⁻¹ night-to-night offsets. To mitigate this observed night-to-night RV offsets induced by the non-stability of the wavelength solutions, a possible solution is to use in combinaison to the Th-Ar spectrum, the much denser spectrum of the FP étalon to stabilise the wavelength solution on the edges of the HARPS CCD (Bauer et al. 2015). This work will be soon implemented in the HARPS DRS (Cersullo et al. 2018) and will be available to the community. Performing a wavelength solution every night simplify a lot the data analysis as this fix the zero velocity point of the instrument every night, and it is therefore not necessary to estimate the drift of the instrument from night-to-night. In Appendix C, we test another solution to mitigate the night-to-night offsets. This solution consists in building a master Th-Ar spectrum from all the Th-Ar calibrations, then applying the wavelength solution of the master to all Th-Ar spectra and measure the drift between them and the master. Doing so, only a unique wavelength solution is used, which solve for the problem of wavelength solution unstability. By using this technique in combination with our new RV extraction approach, we demonstrate that for the 4 years of α Cen B RV data, the standard deviation of the night-to-night systematics can be estimated to 0.4 m s⁻¹. Readers that want more details are invited to have look at Appendix C.

6. Mitigating stellar activity in radial-velocity measurements

Now that we are confident that our new way of extracting the RV information from HARPS spectra gives RVs that are as precise as the ones derived by the HARPS DRS, we can study the RV variation of each individual spectral line.

We want to investigate here the effect of stellar activity on each individual spectral line. Stellar activity, due to the presence of strong magnetic fields in active regions, modifies locally the effective temperature inside spots, and inhibit convection in spots and faculae. We know that each spectral line has a different sensitivity to temperature variation, and therefore, the appearance of a spot on the stellar surface should modify the stellar spectrum. In addition, spectral lines are formed at different depth in the stellar photosphere, and are thus affected differently by convection, as the velocity of convection also varies with depth. In Gray (2009), the author shows that the convective blueshift is of a few hundred of m s⁻¹ for shallow lines, formed deep inside the stellar photosphere, while it is negligible for very deep lines, formed close to the stellar surface. This effect is also well illustrated in the Fig. 3 appearing in Reiners et al. (2016). Therefore, when the convective blueshift is reduced due to the presence of strong magnetic fields in active regions, shallow lines should be significantly affected, while deep lines should not.

Some preliminary works have shown that each spectral line are affected differently by stellar activity. In Davis et al. (2017), the authors show, using PCA on data simulated with SOAP 2.0 (Dumusque et al. 2014), that spectral lines seems to be affected differently by activity. A similar result is also shown in Thompson et al. (2017) and Wise et al. (2018), were the authors demonstrate that stellar activity modifies the equivalent width of spectral lines. However, although variation in spectral line is seen in those analyses, it is not clear if those variations induce a change in RV. Indeed, a spectral line could change in width and/or depth, without changing in asymmetry and therefore without inducing a RV effect.

To analyse the effect of stellar activity on each spectral line, we need RV observations significantly affected by activity, and with a dense enough sampling to detect stellar rotation. We therefore used, like in Thompson et al. (2017), the 2010 data of α Cen B. On those data, we clearly see a ∼5 m s⁻¹ sinusoidal variation due to stellar activity over two rotational periods of the star (Dumusque 2014). A similar signal is also seen in the $\log \left(R_{\mathrm{HK}}^{\prime }\right)$ activity index (Dumusque et al. 2012).

When measuring RVs, the activity signal that we observe is coming from the spectral lines formed inside magnetic regions, faculae and spots, that are formed at the level of the stellar photosphere. The magnetic field that are at the origin of these regions goes up to the chromosphere and create regions that are called plage. The $\log \left(R_{\mathrm{HK}}^{\prime }\right)$ activity index is known to probe plage, and although faculae, spots and plages are related to each other as they originate from the same magnetic fields, it is not necessarily expected that a one-to-one correlation exists between the RV variation induced by stellar activity and the $\log \left(R_{\mathrm{HK}}^{\prime }\right)$ activity index. We therefore decided here to select the RV measured on all the spectral lines as our proxy for the activity signal, and we compared the RV of each individual line relative to it. To mitigate the night-to-night offsets observed in the RVs of spectral lines falling on the edges of the HARPS CCD described in Sect. 5 and in Appendix C, we used the RV data derived with our new RV extraction procedure using a unique wavelength solution. We can see in the right bottom plot of Fig. C.1 that mitigating those night-to-night offsets is extremely important in order to improve the correlation between the RV of each spectral line and the RV measured on all the lines. We finally rejected from this analysis spectral lines falling close to CCD stitchings, to prevent being perturbed by the inaccuracy of the wavelength solution close to those regions (Dumusque et al. 2015).

In Fig. 5, we compare the RVs of three different spectral lines with respect to the activity signal in α Cen B, which is in our case the RVs measured on all the lines (see preceding paragraph). As we can see, the spectral line at 5602.91 Å is strongly correlated with the stellar activity signal, and shows a RV peak-to-peak amplitude of 60 m s⁻¹. On the opposite, we see that the line at 4173.93 Å is anti-correlated with the activity signal. We also observe a lot of lines for which the correlation is null, like for the line at 4337.05 Å. In Fig. 6, we show the median RV error of each spectral line with respect to the R Pearson correlation coefficient between the RVs of those lines and stellar activity. As we can see, we observe 2951 spectral lines that have an absolute correlation weaker than 0.2, in blue, and 489 spectral lines in red that present strong correlation with activity (R > 0.6).

Fig. 5. 2010 RVs of α Cen B as measured on three different spectral lines with wavelength 5602.91 Å (top panel), 4337.05 Å (middle panel) and 4173.93 Å (bottom panel). The median RV error on these spectral lines is 8.7, 8.9 and 10.0 m s−1, respectively. Each panel is divided in two, on the left we show the RV of the line as a function of time, and on the right, the correlation of the RV of this line with the RV derived using all the spectral lines. In 2010, α Cen B is activity and the signal seen in the RV of all the lines is mainly due to activity. Therefore, a strong correlation between the RV of a spectral line with those RVs implies that the line is strongly affect by activity signals.

Now that we observe that the RV of each spectral seems to be affected differently by stellar activity, we want to test if measuring the RV using only affected or non-affected spectral lines can modify the stellar activity signal seen in RV. We therefore calculat the RVs using all the lines available, using the affected red lines shown in Fig. 6, and the non-affected blue lines shown in the same figure. Using the same colour code, we show the result in Fig. 7.

Fig. 6. RV precision on the RVs of each spectral line as a function of the R Pearson correlation coefficient between the RVs of each individual line with the RVs measured on all the lines, as shown in Fig. 5. The blue selection corresponds to spectral line for which the absolute R correlation coefficient is below 0.2, and therefore to lines that does not seem to be strongly affected by stellar activity. The red select includes lines for which the correlation is stronger than 0.6, and therefore lines that are sensitive to activity.

When combining the RVs of the 489 spectral lines strongly affected by stellar activity, we obtain RVs that present a standard deviation that is twice higher than the one measured on the RVs of all the lines. By using the RV error of each of these lines, which is shown in Fig. 6, we can compute the RV weight of those 489 spectral lines with respect to the RV weight of all the spectral lines, which by definition is 100%. Doing so, we find that the RV weight of those 489 lines is 16.5%, which implies that the photon noise will be $1/\sqrt{0.165}=2.5$ times larger compared to the photon noise obtained when measuring the RV on all spectral lines. For α Cen B, the RV photon noise is so low that this increase by a factor 2.5 will not have a significant impact, however it could be the case for fainter stars.

When combining the RVs of the 2951 spectral lines that show a weak correlation with activity, we obtain the blue RVs in Fig. 7. The standard deviation of those RVs is 1.6 times smaller than the standard deviation measured on the RVs derived using all the spectral lines. Therefore, the activity signal is significantly mitigated when using this specific selection of spectral line. Those 2951 lines contain 26.6% of the total RV information, implying that the photon noise on the RVs obtained from this line selection is $1/\sqrt{0.266}=1.9$ times larger than the photon noise obtained when measuring the RV on all spectral lines. By measuring the RVs only on spectral lines not significantly affected by activity, we are able to mitigate the red noise induced by stellar activity by a factor of 1.6 while increasing the white photon-noise by a factor of 1.9. Mitigating the red noise induced by stellar activity has more weight than increasing the white photon noise when searching for exoplanet signal in RV data.

We are confident that deriving the RVs using sub-samples of spectral lines is the key to mitigate the stellar activity signal. When looking at the RV derived using only unaffected lines, it seems that we still see some signal from stellar activity. What we could do is looking at the correlation between the RV of the affected and unaffected lines, and removing it to mitigate even better the activity signal in the unaffected RVs. We however have to be careful, because if a planet is present in the data, its signal will be the same in the RV of the affected or unaffected lines. By removing the correlation, we could modify the planetary signal in the RV residuals and thus obtain incorrect orbital parameters. In conclusion, we have to find what is the best technique to correct for stellar activity using the RV measured on sub-samples of spectral lines, and this is what we will investigate in a forthcoming paper (Cretignier et al., in prep.).

Fig. 7. Comparison of the 2010 RVs of α Cen B as derived when using all the spectral line (in black), when using lines for which the correlation between their RVs and the RV of all the lines is strong (in red) and when using lines for which the correlation is weak (in blue). The blue and red selections of lines can be see in Fig. 6. We highlight as well, in the title of each subplot, the number of lines used in each selection, and by how much those lines contribute to the RVs derived using all the spectral lines available. We clearly see that by selecting lines that are strongly correlated or not with the activity signal, we are able to boost the activity signal seen in the RVs by a factor of 2, or mitigate it by a factor of 1.6.

7. Discussion

In this paper, we present a new approach to derive RVs from high-resolution spectra. First, we measure the RV of each individual spectral line and then we combine the RV information of all the spectral lines together to obtain a precise stellar RV.

In Sect. 4 we compare the RVs derived by our new RV extraction procedure, with the RVs obtained from the HARPS DRS, which is the gold standard for precise RV estimate from G and K dwarfs HARPS spectra. For the three stars studied, τ Ceti, α Cen B and HD 10180, we find very similar results (see Figs. 1, 3 and 4). In the case of τ Ceti (see Sect. 4.1), although the RVs are extremely similar when looking at the RV difference between the two reduction and at the distribution of the derived RVs, it seems that ou new RV extraction procedure is slightly more sensitive to systematics that create spurious RV signals between 100 and 1000 days in period. Those signals are however not significant when considering a false alarm probably of 1%. We still try to investigate the origin of those extra signals in Appendix B, were we show that they might be due to the fact that our reduction is more sensitive to systematics induced by spectral line crossing CCD stitchings, and by the effect of tellurics. For HD 10180 (see Sect. 4.3), besides checking that the RVs from both reductions are extremely similar, we also compared the orbital parameters derived for the planets present in the data. Using the two different RV reductions, we found exactly the same orbital parameters considering their 1-σ uncertainty. Therefore, the RVs derived with our new RV extraction procedure gives RVs that are as good as the HARPS DRS to search for exoplanets.

Using our new RV extraction procedure, we also estimated the drift of the spectrograph overnight, by measuring the RV as a function of time for each of the Th-Ar or FP étalon emission spectral lines that appear on the simultaneous reference spectra that are taken contemporaneously with stellar observations. As we can see in Fig. A.1, our new RV extraction procedure gives also extremely similar results than the HARPS DRS when estimating the drift of the instrument overnight.

Looking at the RV of each individual spectral line, we saw in Sect. 5 and in Appendix C that the RVs of lines falling on the edges of the HARPS detector can be strongly affected by night-to-night offsets. We found that these RV offsets are induced by using a different wavelength solution every night.

To correct for these night-to-night RV offsets, we choose in this paper to build a master Th-Ar spectrum and to apply to each Th-Ar calibrations the unique wavelength solution of this Th-Ar master. This prevent us of using non-stable wavelength solutions every night, however this requires to correct for the drift of the instrument over time. Although this approach was significantly complicating the data reduction process (see Appendix C), we were able to demonstrate that this solution can mitigate the night-to-night RV offsets seen on spectral lines falling on the edges of each spectral order. By deriving the RV of α Cen B from 2008 to 2012 using a unique wavelength solution, we could estimate that the standard deviation of the night-to-night instrumental systematics is 0.4 m s⁻¹. The obtained RVs are therefore less sensitive to the night-to-night RV offsets, and the activity signal of the star, seen at the rotational period, is better characterised. We however discuss in Sect. 5 and Appendix D that due to a slow modification of the spectrograph focus over time, the RVs derived using a unique wavelength solution will be affected by a long-term drift, that can however be correct for.

Although we could show that our new RV extraction procedure using a unique wavelength solution reduces the night-to-night RV variations, the final gain in terms of RV precision is rather small. Therefore, what is the gain of doing this complicated reduction, while a simple CCF approach gives extremely similar results.

First, the main advantage, and this is the most relevant astrophysical phenomena that our new RV extraction demonstrate, is the observation that some lines are much more affected by stellar activity effects than others (see Sect. 6). This is expected from stellar physics, however it was never demonstrated that using current spectrographs used for exoplanet characterisation, we had the precision to measure solar-like stellar activity in the RVs of individual spectral lines. By measuring the RV of α Cen B using only lines strongly affected or unaffected by stellar activity, we were able to boost the activity signal in RVs by a factor of 2, or mitigate it by a factor of 1.6, respectively. Deriving the RVs on a subset of spectral lines has for consequence the increase in photon noise of the final RVs. Using the line selection that doubles or mitigate by a factor of 1.6 stellar activity worsen the RV precision by a factor of 2.5 or 1.9, respectively. On α Cen B, this does not have any consequences as the photon noise for this target is extremely low, however, it might be a problem for fainter targets. This is something that we will investigate soon, as we believe that for planet detection and characterisation, mitigating the red noise induced by stellar activity has much more weight than increasing the white photon noise of the derived RVs.

In addition, our new RV extraction approach can be used to probe how the RV varies locally on the detector, and not only from orders to orders, what the HARPS DRS and other data reduction package do, but also within orders. From Figs. D.1 and D.2 it is clear that the drift measured on the instrument over time is not the same at every location on the HARPS CCD, and we see significant differences between the blue and red detectors, but also between the left and right sides of the orders. Our new RV extraction technique is therefore interesting to probe and understand instrumental systematics.

Moreover, it is important to note that once the RV of each line is computed, it takes nearly no time to measure the RV for any given selection of line, as this process only requires a weighted average. This is of course extremely interesting if we want to compare the RV derived from a huge number of different line selections.

8. Conclusion

In this paper, we present a new approach to derive precise RVs from HARPS spectra. This new method tries to use the wealth of information contained in a high-resolution spectrum rather than averaging all the information into the few parameters of the CCF. We demonstrate that this new RV extraction procedure can be used to probe stellar spectral lines that are sensitive or not to activity, and that by using a smart selection of them, we can either boost or mitigate the stellar activity signal observed in RV. This is of course very interesting in the light of new spectrographs like ESPRESSO that will reach a precision of 0.1 m s⁻¹ and for which stellar activity signal will dominate the data. We will therefore investigate in more details the effect of activity on each spectral line and try to understand the underlying physics in a forthcoming paper (Cretignier et al., in prep.). This forthcoming work will present in detail the lines that are strongly affected or unaffected by stellar activity, with their physical properties. We hope that this will motivate other people in the community to join us in this effort for solving the problem of stellar activity in RV measurements, which would be beneficial for the entire exoplanet community.

This new RV extraction procedure can also be used to study the impact of tellurics on the RVs (see Sect. B.2 and Cunha et al. 2014), but also to understand better instrumental systematics (see for example Figs. D.1 and D.2).

As a ending note, we want to emphasise that our new RV extraction procedure demonstrated that deriving the RV on a spectral line basis brings a huge amount of information, mainly to understand better and correct for stellar activity, and thus improving the methods for extracting locally the RV in a stellar spectra should be investigated further. The results shown here are not specific to HARPS, and can be obtained for the data of any stabilised high-resolution spectrographs.

Acknowledgments

XD would like to thank F. Pepe, C. Lovis, F. Bouchy and D. Segransan for their help in understanding all the details of the HARPS spectrograph and its data reduction software. XD is grateful to The Branco Weiss Fellowship–Society in Science for its financial support. This work has been carried out in the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). This research has made use of the DACE platform developed in the framework of PlanetS (https://dace.unige.ch).

References

Appendix A: Correcting for the HARPS instrumental drift during the night

Before the beginning of each night, the HARPS science fibre (fibre A) is illuminated with a Th-Ar lamp to perform a wavelength solution and thus fix the zero point of the instrument. At the same time as the Th-Ar spectrum is recorded on the science fibre, the spectrum of another lamp, either Th-Ar or FP étalon, is recorded on the HARPS reference fibre (fibre B). During the night, this second lamp continues to illuminate the reference fibre, while stars are observed on the science fibre. This allows to probe the instrumental drift within a night by measuring the drift between this simultaneous reference spectrum and the calibration made before the night. In the HARPS DRS, this drift is measured using the technique described in Sect. 3.4 comparing each order of the Th-Ar or FP simultaneous spectra with the respective order in the Th-Ar or FP calibration. At the end, a weighted average is performed on all orders to get a measure of the instrumental RV drift for each observation. Here, for consistency with the technique used to derive the RV on each individual line, we build a Th-Ar and/or FP master spectrum from all the simultaneous spectra taken during one night and then measure, on each emission line of the Th-Ar or FP simultaneous spectra, the RV drift relative to the Th-Ar or FP master, respectively. Once this is done, we can either perform a weighted average on all the spectral lines to get a single RV drift for each observation, like it is done in the HARPS DRS, or we can measure locally the drift on the HARPS detector.

In Fig. A.1, we compare the drifts estimated by the HARPS DRS and the drifts derived by our new RV extraction procedure using all the simultaneous spectra obtained for the τ Ceti RV data set that we analyse in Sect. 4.1. As we can see, our new RV extraction procedure gives extremely similar results than the HARPS DRS and for consistency we use those newly derived drifts to correct for the HARPS systematics when we estimate the RV of each individual spectral line.

Appendix B: Signals between 100 and 1000 days in τ Ceti RVs

As we can see in Fig. 1 of Sect. 4.1, the RVs derived by our new RV extraction procedure seems to be more affected by spurious RV signals with periods between 100 and 1000 days. Our first guess was that our reduction was more sensitive to the systematics created by spectral lines falling close to CCD stitchings, as it was shown by Dumusque et al. (2015) that those lines are responsible for inducing spurious RV signal with periods of one-year and its different harmonics. We therefore removed 1234 spectral lines that were too close to CCD stitchings and recomputed, using our new RV extraction procedure and the HARPS DRS, the RVs for τ Ceti. The results are shown in Fig. B.1. We see that some of the strong signals between 100 and 1000 days disappear in this new analysis, however, not all of them. The behaviour of spectral lines close to CCD stitchings can therefore not explain all the observed systematics.

Fig. A.1. Top panel: comparison of the drift of the instrument measured on the simultaneous reference fibre illuminated either by the Th-Ar or the FP étalon lamp and estimated by our new RV extraction procedure (in red) and by the HARPS DRS (in green). Bottom panel: histogram of the drifts measured. The legend shows the standard deviation obtained for both drift data sets.

Fig. B.1. Fig. 1 but for the RVs of τ Ceti derived without considering spectral lines that falls close to CCD stitchings.

Fig. B.2. RVs of τ Ceti derived only using the spectral lines with wavelength longer than 5000 Å (left panel), shorter than 5000 Å (middle panel) and shorter than 5000 Å plus a quadratic drift fitted to the residual RVs to remove low frequencies (right panel). This quadratic drift corresponds to a small variation of the RVs of less than 1.5 m s−1 over 11 years. It is clear that spectral lines with wavelengths longer than 5000 Å induce significant signals at periods between 100 and 1000 days. This is certainly due to telluric and micro-telluric lines affecting the red part of the spectrum.

Another phenomenon that can create signals with a one-year period is the one caused by tellurics. Tellurics do not move on the HARPS CCD, while stellar spectra are red and blueshifted over a year due to the revolution of Earth around the Sun, which modifies the RV of the Earth in the direction of the observed star. Tellurics can therefore induce signals with periods of a year and its different harmonics, depending on the sampling of the data. To investigate if tellurics could be at the origin of the yearly signal observed, we compared the RVs of all the lines with a wavelength smaller than 5000 Å with the RVs of all the lines at longer wavelengths. We selected a cutoff at 5000 Å here as there is no significant tellurics below this wavelength. In Fig. B.2, we can see that for the lines redder than 5000 Å, significant signals appear between 100 and 1000 days, while for the lines bluer than 5000 Å, no significant signal is observed, except a long-term trend. It seems therefore that the extra signals that we see in our new RV extraction procedure RVs at periods around 100 and 1000 days compared to the HARPS DRS RVs is induced by spectral lines redder than 5000 Å. Because this is where telluric contamination is strong, it is likely that those extra signals are induced by tellurics. This is however something that we will investigate with more precision in the future.

Appendix C: Correcting night-to-night RV offsets due to variation in wavelength solutions

As explained in Sect. 5, a wavelength solution is derived every night to fix the zero velocity point of the HARPS spectrograph. To derive a wavelength solution on a Th-Ar calibration, a third order polynomial is fitted on the Th-Ar emission lines present in each order. Because the flux is low on the edges of the HARPS orders due to the blaze, the polynomial will not be constrained at those locations and the wavelength solution can vary strongly from night-to-night.

To show the non-stability of the wavelength solution on the edges of the HARPS orders, we plot in the left panel of Fig. C.1 the differences in RV between consecutive wavelength solutions for HARPS spectral order 38 obtained in 2010. For reasons that will become clear in Sect. 6, we analysed only the wavelength solutions used for the RV data of α Cen B in 2010. As we can see, for pixel values larger than 3500, we observe two categories of wavelength solutions: one category that presents negative RVs, and another one that exhibits positive RVs. The difference in RV between those two categories of wavelength solution for pixels values larger than 3500 is ∼15 m s⁻¹. With our new RV extraction procedure, we can study the behaviour of a spectral line falling in this region of the detector. We therefore analysed the variation of the 4999.5 Å Ti spectral line, which is localised at pixel 3800 in order 38 (red vertical dashed line in the left panel of Fig. C.1). In the top right plot of Fig. C.1, we show in green the RVs of this spectral line as a function of time as measured with our new RV extraction procedure. Like for the wavelength solutions, we observe that the RVs can also be classified in two categories: one with positives RVs and another one with negative RVs. The difference between those two categories being also ∼15 m s⁻¹. Our new RV extraction procedure therefore allows us to see directly the impact of the non-stability of the wavelength solution on the RVs derived for the 4999.5 Å Ti spectral line.

Fig. C.1. Left panel: wavelength solutions used for the 2010 data of α Cen B, for spectral order 38. As we can see, those wavelength solutions are not stable on the edges of order 38. Top right panel: RV measured on the 4999.51 Å Ti line that falls on pixel 3800 of order 38 (red vertical dashed line in the left panel). In green, we show the RVs as measured when using the wavelength solutions shown on the left panel. In red we show the RVs derived when using a unique wavelength solution. As we can see, the non-stability of the wavelength solutions on the right edge of spectral order 38 induces night-to-night RV offsets of ∼15 m s−1, which is mitigated when using a unique wavelength solution. Bottom right panel: correlation between the RV of the 4999.51 Å Ti line and the RV derived using all the spectral lines. As we can see, the correlation is significantly stronger when we use the RVs derived using a unique wavelength solution.

To mitigate this observed night-to-night RV offsets induced by the non-stability of the wavelength solutions, a possible solution is to use in combinaison to the Th-Ar spectrum, the much denser spectrum of the FP étalon to stabilise the wavelength solution on the edges of the HARPS CCD (Bauer et al. 2015). This work will be soon implemented in the HARPS DRS (Cersullo et al. 2018) and will be available to the community. Performing a wavelength solution every night enables to fix the zero velocity point of the instrument, which simplify greatly the data reduction as it is not necessary to estimate the drift of the instrument from night-to-night. If we want to avoid using a wavelength solution every night, we have to measure the drift of the instrument from night-to-night. This can be done by building a master Th-Ar spectrum from all the Th-Ar calibrations, then applying the wavelength solution of the master to all Th-Ar spectra and measure the drift between them and the master. This is the method that we explore here and for simplicity, we will refer to this technique in the rest of the paper as the one using a unique wavelength solution.

This approach complicates somehow the data analysis. First, within the 15 years of the HARPS spectrograph, the Th-Ar lamp used for deriving wavelength solutions was changed five times, and thus we need to create a Th-Ar master spectrum for each of these lamps. In addition, several interventions on the instrument modified the way the Th-Ar spectrum is recorded on the detector. On JD = 24571745 (1st of June 2015), the fibres to guide the light from the ESO 3.6 m telescope to HARPS where changed from circular to octagonal, to mitigate the RV effect induced by telescope guiding errors. This drastically changed the point spread function (PSF) of the spectrograph, which induced a RV offset of ∼2 km s⁻¹ (see Fig. D.1). Therefore, a new Th-Ar master spectrum needed to be created for analysing the data taken after the fibre change, even though the Th-Ar lamp did not change. On JD = 2453399 (19th of January 2005), an intervention on the instrument modified the relative drift between the blue and red CCD of HARPS, which can be seen in Fig. D.2. Besides a general offset and a relative difference between the blue and red detectors of HARPS, we can also see in Figs. D.1 and D.2 that the RV offset induced by those interventions is larger on the left side than the right side of the CCD. All those local differences prevent us of measuring a single RV drift for the entire detector to correct for the drift of the spectrograph over time. To solve for this problem, we divided the HARPS CCD in 24 regions, 16 corresponding to the blue CCD (pixels 0–1024, 1024–2048, 2048–3072, 3072–4096 and orders 0 to 12, 12–24, 24–36 and 36–47) and 8 corresponding to the red CCD (same separation in pixels but for orders 47–60 and 60–72). We then measured the drift of the instrument in each of these regions by averaging the RV of all the lines falling in these regions. This allows to correct locally for the drift of the spectrograph over time.

Fig. C.2. Ratio between the standard deviations of the RVs of each spectral line measured using one wavelength solution and using a wavelength solution every night. As we can see, on the edges of the HARPS CCD, which corresponds to pixels smaller than 1000 and larger than 3000, our new RV extraction procedure using only one wavelength solution give RVs for which the standard deviation is smaller than when deriving the RVs using a wavelength solution every night.

Fig. C.3. Same as Fig. 1 but for the RVs of α Cen B using only one wavelength solution.

By measuring the drift of the instrument using a unique wavelength solution, we strongly reduce the night-to-night RV offsets observed on spectral lines falling on the edge of the HARPS CCD. In the top right plot of Fig. C.1, we show in red the RVs of the 4999.5 Å Ti spectral line derived with this new approach. As we can see, the night-to-night RV offsets are strongly mitigated, which is shown by the reduction of the RV standard deviation from 8 to 3.4 m s⁻¹. On the bottom right plot, we see how the correlation between the RVs of the Ti line and the RVs of all the lines become stronger when we use a unique wavelength solution. This will be very important for the analysis performed in Sect. 6.

To illustrate that using a unique wavelength solution helps in reducing the night-to-night RV offsets of all the lines falling on the edges of the HARPS CCD, we calculated the ratio between the standard deviations of the RVs measured using a unique wavelength solution and the RVs derived using a wavelength solution every night. The result for all the lines present in the 2010 spectra of α Cen B is shown in Fig. C.2. As we can see, this ratio is on average always smaller than one, and the deviation from one becomes more and more important towards the edges of the HARPS CCD. This implies that by measuring the RV drift of the instrument using a unique wavelength solution, we can mitigate significantly the night-to-night RV offsets induced by the non-stability of the wavelength solution computed every night.

Fig. C.4. Histogram of the absolute night-to-night RV offsets measured on the 2008–2012 data of α Cen B. In green, we represent the night-to-night offsets as measured in the RVs from the HARPS DRS, and in red we show the same offsets obtained from the RVs derived using a unique wavelength solution. As we can see, the RV standard distribution of the red histogram is smaller, proving that the night-to-night offsets are reduced by our analysis.

In this section, we show that measuring the instrument drift over time using a unique wavelength solution gives RVs that are less affected by night-to-night systematics compared to using a wavelength solution every night. In Fig. C.3, we compare the RVs of α Cen B estimated from the HARPS DRS and derived with our new RV extraction procedure using a unique wavelength solution. As we can see, those two sets of RVs are extremely similar, with a standard deviation of the RV differences of 0.59 m s⁻¹. Looking at the periodogram of those RVs, we see that the activity signal seen at the stellar rotation period, near 40 days, is stronger in our new RV extraction procedure. This comes from the fact that by mitigating the night-to-night RV offsets using a unique wavelength solution, the activity signal becomes better characterised. The mitigation of this systematics can also be seen in Fig. C.4, were we plot the histogram of the night-to-night RV offsets as measured on the RV data of α Cen B. As we can see, the standard deviation of the absolute night-to-night RV offsets is smaller when using a unique wavelength solution. Noises adds up quadratically, therefore this analysis shows that for the RV data of α Cen B, the noise induced by the non-stability of the wavelength solutions is 0.4 m s⁻¹ $(=\sqrt{{1.56}^2-{1.51}^2})$.

Fig. D.1. Drift of HARPS as a function of time measured on the Th-Ar spectra taken each night to derive the wavelength solution of the detector. Each subplot represents the drift in RVs measured on the Th-Ar emission lines falling in different regions of the HARPS detector. The HARPS detector is cut into 24 regions, 16 corresponding to the blue CCD shown with a blue background (pixel 0–1024, 1024–2048, 2048–3072, 3072–4096 and orders 0–12, 12–24, 24–36 and 36–47) and 8 corresponding to the red CCD shown with a red background (same separation in pixel but for order 47–60 and 60–72). The main feature that we see here is the huge RV offset on JD = 2457174 (1st of June 2015), which is induced by the change of the HARPS fibres from circular to octagonal.

In the top plot of Fig. C.3, we corrected the RVs from the long-term trend induced by the presence of α Cen A using a second order polynomial. The final RVs using a unique wavelength solution are extremely similar to the HARPS DRS RVs after this drift correction, however, the fitted drift is not the same. As this drift is induced by α Cen A, it should be the same. Therefore, the only explanation for this difference is that deriving the RVs using a unique wavelength solution introduce a long-term drift in the RVs. This extra drift is due to the change in focus of the HARPS spectrograph with time. We are able to measure this effect in Appendix D, and correct for it.

As a note of caution, if a FP étalon is available at a RV facility, we strongly encourage the use of a wavelength solution every night if the FP spectrum can be used in combination with the Th-AR spectrum to derive a stable wavelength solution (Bauer et al. 2015; Cersullo et al. 2018). We demonstrated here that using a single wavelength solution works as well, however it complicates quite significantly the data analysis process.

Appendix D: HARPS systematics

By measuring with our new RV extraction procedure the RVs on each emission spectral line in the Th-Ar spectra used to perform the wavelength solution every night, we can study the drift of the HARPS spectrograph over time. In Fig. D.1, we show the drift observed when the fibre from HARPS were changed from circular to octagonal on the 1st of June 2015. As we can see, this intervention on the instrument induced a RV offset of more than 2 km s⁻¹. In addition, by looking carefully, we can also see that the offset is not the same on the right side and on the left side of the detector, with a difference of about 50–100 m s⁻¹.

Fig. D.2. Same as Fig. D.1 but here we zoom on the instrument offset happening on JD = 2453390 (19th of January 2005).

In Fig. D.2, we show the drift of the instrument induced by an intervention on HARPS on the 19th of January 2005 (JD = 2453390). In this case, we see a different RV offset for the blue and the red CCD, ∼100 and ∼10 m s⁻¹, respectively. In this case, we also see a difference in RV offset between the left and right side of the HARPS CCD.

Fig. D.3. Same configuration as Fig. D.1. This time however, we show the drift between the Th-Ar spectra with their own wavelength solution with respect to Th-Ar master spectra. We note that we have a master Th-Ar spectrum for each of the 5 lamps used over the 15 years of HARPS, plus an additional one to compensate for the systematics induced by the change of the HARPS fibres. In this plot, we corrected for the offset that exists between each consecutive Th-Ar master spectra. The drift observed over time is due to a slow variation of the focus of the instrument, which changes the PSF and therefore the shape of the spectral lines over time in the Th-Ar spectra recorded on the CCD. The master Th-Ar spectra is by construction an average of all the Th-Ar spectra, and therefore do not account for the change in the focus of HARPS.

To correct for the night-to-night RV offsets due to the non-stability of the wavelength solution on the edges of HARPS orders, we first built a master Th-Ar spectrum from all the Th-Ar calibration spectra. Then, we applied for all Th-Ar calibrations the wavelength solution of the master Th-Ar spectrum and measured the drift between each of these calibrations and the master. We note that because we used five Th-Ar calibrations lamps over the 15 years of HARPS, we needed to build a master Th-Ar spectrum for each of them. In addition, an extra master was also considered to account for the drastic change induced by the modification of the HARPS fibres. Besides changing the calibration lamp that induces a RV offset but which is easy to correct for, this method requires that the spectrum of the Th-Ar lamp recorded on the CCD does not change over time. On HARPS, we know that this is not true as that the focus of the instrument slowly varies with time. This implies that the PSF of the instruments changes, and therefore the shape of the lines in the Th-Ar spectrum recorded on the detector changes as well. To measure this effect, we estimated the drift between each Th-Ar calibration with their own wavelength solution and the master Th-Ar spectra used. The result of this analysis is shown in Fig. D.3. If the shape of the Th-Ar spectral lines varies with time, this analysis should reveal a long-term drift. This is exactly what is seen, with a drift of ∼50 m s⁻¹ over 15 years. This proves that the focus of the instrument is slightly changing over time. From this result, we can conclude that using a unique wavelength solution will induce a drift in the derived RVs. This drift can be corrected for by using the results shown in Fig. D.3, or simply by fitting a polynomial as we did for the RVs of α Cen B in Sect. 5.

All Tables

All Figures

Fig. 1. Top panel: comparison of the RVs of τ Ceti derived using our new RV extraction procedure (red) with the RVs estimated from the HARPS DRS (green). The red vertical dashed line corresponds to the change from circular to octagonal fibres that was performed on HARPS on the 1st of June 2015 (JD = 2457174). Middle panel: difference between the RVs derived using our new RV extraction procedure and the RVs from the HARPS DRS. We show in the legend the standard deviation of the RV difference for the measurements preceding the fibre change. Bottom panel: periodograms of the RVs derived using our new RV extraction procedure (red) and the RVs derived by the HARPS DRS (green). To calculate those periodograms, we rejected the data after the fibre change to prevent low frequency signals induced by the observed offset. The two periodograms have been normalised in such a way that their FAP of 1% coincide.

In the text

	Fig. 2. Histogram of the RVs of τ Ceti, HD 128621 and HD 10180 derived using our new RV extraction procedure (red) and the RVs derived using the HARPS DRS (green). The legend shows the standard deviation obtained for both data sets.
In the text

	Fig. 3. Same as Fig. 1 but for the RVs of α Cen B.
In the text

	Fig. 4. Same as Fig. 1 but for the RVs of HD 10180.
In the text

Fig. 5. 2010 RVs of α Cen B as measured on three different spectral lines with wavelength 5602.91 Å (top panel), 4337.05 Å (middle panel) and 4173.93 Å (bottom panel). The median RV error on these spectral lines is 8.7, 8.9 and 10.0 m s−1, respectively. Each panel is divided in two, on the left we show the RV of the line as a function of time, and on the right, the correlation of the RV of this line with the RV derived using all the spectral lines. In 2010, α Cen B is activity and the signal seen in the RV of all the lines is mainly due to activity. Therefore, a strong correlation between the RV of a spectral line with those RVs implies that the line is strongly affect by activity signals.

In the text

Fig. 6. RV precision on the RVs of each spectral line as a function of the R Pearson correlation coefficient between the RVs of each individual line with the RVs measured on all the lines, as shown in Fig. 5. The blue selection corresponds to spectral line for which the absolute R correlation coefficient is below 0.2, and therefore to lines that does not seem to be strongly affected by stellar activity. The red select includes lines for which the correlation is stronger than 0.6, and therefore lines that are sensitive to activity.

In the text

Fig. 7. Comparison of the 2010 RVs of α Cen B as derived when using all the spectral line (in black), when using lines for which the correlation between their RVs and the RV of all the lines is strong (in red) and when using lines for which the correlation is weak (in blue). The blue and red selections of lines can be see in Fig. 6. We highlight as well, in the title of each subplot, the number of lines used in each selection, and by how much those lines contribute to the RVs derived using all the spectral lines available. We clearly see that by selecting lines that are strongly correlated or not with the activity signal, we are able to boost the activity signal seen in the RVs by a factor of 2, or mitigate it by a factor of 1.6.

In the text

	Fig. A.1. Top panel: comparison of the drift of the instrument measured on the simultaneous reference fibre illuminated either by the Th-Ar or the FP étalon lamp and estimated by our new RV extraction procedure (in red) and by the HARPS DRS (in green). Bottom panel: histogram of the drifts measured. The legend shows the standard deviation obtained for both drift data sets.
In the text

	Fig. B.1. Fig. 1 but for the RVs of τ Ceti derived without considering spectral lines that falls close to CCD stitchings.
In the text

Fig. B.2. RVs of τ Ceti derived only using the spectral lines with wavelength longer than 5000 Å (left panel), shorter than 5000 Å (middle panel) and shorter than 5000 Å plus a quadratic drift fitted to the residual RVs to remove low frequencies (right panel). This quadratic drift corresponds to a small variation of the RVs of less than 1.5 m s−1 over 11 years. It is clear that spectral lines with wavelengths longer than 5000 Å induce significant signals at periods between 100 and 1000 days. This is certainly due to telluric and micro-telluric lines affecting the red part of the spectrum.

In the text

Fig. C.1. Left panel: wavelength solutions used for the 2010 data of α Cen B, for spectral order 38. As we can see, those wavelength solutions are not stable on the edges of order 38. Top right panel: RV measured on the 4999.51 Å Ti line that falls on pixel 3800 of order 38 (red vertical dashed line in the left panel). In green, we show the RVs as measured when using the wavelength solutions shown on the left panel. In red we show the RVs derived when using a unique wavelength solution. As we can see, the non-stability of the wavelength solutions on the right edge of spectral order 38 induces night-to-night RV offsets of ∼15 m s−1, which is mitigated when using a unique wavelength solution. Bottom right panel: correlation between the RV of the 4999.51 Å Ti line and the RV derived using all the spectral lines. As we can see, the correlation is significantly stronger when we use the RVs derived using a unique wavelength solution.

In the text

Fig. C.2. Ratio between the standard deviations of the RVs of each spectral line measured using one wavelength solution and using a wavelength solution every night. As we can see, on the edges of the HARPS CCD, which corresponds to pixels smaller than 1000 and larger than 3000, our new RV extraction procedure using only one wavelength solution give RVs for which the standard deviation is smaller than when deriving the RVs using a wavelength solution every night.

In the text

	Fig. C.3. Same as Fig. 1 but for the RVs of α Cen B using only one wavelength solution.
In the text

Fig. C.4. Histogram of the absolute night-to-night RV offsets measured on the 2008–2012 data of α Cen B. In green, we represent the night-to-night offsets as measured in the RVs from the HARPS DRS, and in red we show the same offsets obtained from the RVs derived using a unique wavelength solution. As we can see, the RV standard distribution of the red histogram is smaller, proving that the night-to-night offsets are reduced by our analysis.

In the text

Fig. D.1. Drift of HARPS as a function of time measured on the Th-Ar spectra taken each night to derive the wavelength solution of the detector. Each subplot represents the drift in RVs measured on the Th-Ar emission lines falling in different regions of the HARPS detector. The HARPS detector is cut into 24 regions, 16 corresponding to the blue CCD shown with a blue background (pixel 0–1024, 1024–2048, 2048–3072, 3072–4096 and orders 0–12, 12–24, 24–36 and 36–47) and 8 corresponding to the red CCD shown with a red background (same separation in pixel but for order 47–60 and 60–72). The main feature that we see here is the huge RV offset on JD = 2457174 (1st of June 2015), which is induced by the change of the HARPS fibres from circular to octagonal.

In the text

	Fig. D.2. Same as Fig. D.1 but here we zoom on the instrument offset happening on JD = 2453390 (19th of January 2005).
In the text

Fig. D.3. Same configuration as Fig. D.1. This time however, we show the drift between the Th-Ar spectra with their own wavelength solution with respect to Th-Ar master spectra. We note that we have a master Th-Ar spectrum for each of the 5 lamps used over the 15 years of HARPS, plus an additional one to compensate for the systematics induced by the change of the HARPS fibres. In this plot, we corrected for the offset that exists between each consecutive Th-Ar master spectra. The drift observed over time is due to a slow variation of the focus of the instrument, which changes the PSF and therefore the shape of the spectral lines over time in the Th-Ar spectra recorded on the CCD. The master Th-Ar spectra is by construction an average of all the Th-Ar spectra, and therefore do not account for the change in the focus of HARPS.

In the text