New wavelength calibration of the HARPS spectrograph

TL;DR

This paper presents an improved wavelength calibration method for the HARPS spectrograph, reducing systematic errors and enhancing its ability to detect low-amplitude exoplanet signals by correcting CCD stitching anomalies and using a new thorium line list.

Contribution

It introduces a new calibration approach that accounts for CCD pixel anomalies and employs a higher-accuracy thorium line list, improving radial velocity precision for HARPS data.

Findings

Reduced radial velocity dispersion in HARPS data.

Attenuated one-year period signals in stellar RV measurements.

Enhanced detection sensitivity for low-amplitude exoplanet signals.

Abstract

(abridged) Even if the HARPS spectrograph has been operational for more than 15 years and it provides among the most precise Doppler measurements, improvements are still possible. One known problem, for instance, is the non-fully regular block-stitching of the CCDs, which introduces, in some cases, one-year period parasitic signals in the measured radial velocity. The aim is to improve the wavelength calibration of HARPS to push further its planet-detection capabilities. The properties of the CCD stitching-induced pixel-size anomalies are determined with LED flat-field frames, and then a physical, gap-corrected map of the CCDs is used for the fitting model of the spectral orders. We also use a new thorium line list, based on much higher-accuracy measurements than the one used up to now. We derive new wavelength solutions for the 15 years of HARPS data, both before and after the 2015 fibre upgrade. We demonstrate that we correct the gap anomalies by computing the wavelength solutions of laser frequency comb exposures, both with and without taking the gap correction into account. By comparing the rms of the most stable stars of the HARPS sample, we show that we globally decrease the radial velocity dispersion of the data, especially for the data acquired after the change of fibres. Finally, the comparative analysis of several individual systems shows that we manage to attenuate the periodogram power at one year in most cases. The analysis of the RVs derived from individual stellar lines also shows that we correct the stitching-induced RV variation. This improved calibration of the HARPS spectrograph allows to go deeper in the search for low-amplitude radial-velocity signals. It will be further improved by combining the thorium calibration spectra with laser frequency comb and Fabry-Perot calibration spectra, and not only for HARPS but notably also for HARPS-N and ESPRESSO.

Tables

Table 1. Table 1: Estimated gap sizes. L nom = 15   µ ​ m subscript 𝐿 nom times 15 micrometer L_{\mathrm{nom}}=$15\text{\,}\mathrm{\SIUnitSymbolMicro m}$ .

block	Estimated gap sizes [$\times {10}^{-3}{L}_{\mathrm{nom}}$]
group	1	2	3	4
1	-9.92	-19.67	-29.72	-26.34
	-18.28	-30.71	+3.33	-16.12
2	-5.69	+13.48	+5.20	+8.97
	+1.39	+12.28	-14.36	+2.96
3	-18.30	-9.14	-11.02	-1.57
	-19.62	-16.73	-4.16	- 6.99
4	+0.88	-12.49	+5.02	-16.08
	+0.28	-17.63	-15.62	-21.25
5	+0.26	-1.92	-9.77	+3.45
	-8.46	-5.86	-4.52	+1.90
6	-23.84	-17.08	-27.46	-13.80
	-29.01	-24.26	-34.54	-32.18
7	-15.79	+6.88	+5.04	+15.27
	-5.77	-6.46	-8.16	+6.38

Table 2. Table 2: Comparison of the standard (std) and the new versions of the DRS. The numbers of points are the total and, below, the number before+after the 2015 upgrade. All these points are night-averaged. The rms columns provide the post-fit RV dispersions, both before and after the fibre upgrade. All the rms are in m   s − 1 times meter second 1 \mathrm{m}\text{\,}{\mathrm{s}}^{-1} . The timespans Δ ​ t Δ 𝑡 \Delta t are in days. Arrows mean “std value” → → \rightarrow “new value”. “ n 𝑛 n  pl.” means n 𝑛 n planetary (Keplerian) signals; “MC” means “magnetic cycle”, also fitted by a Keplerian signal unless stated otherwise. References: D16 = Díaz et al. ( 2016 ) ; F17a = Feng et al. ( 2017a ) ; F17b = Feng et al. ( 2017b ) ; L06 = Lovis et al. ( 2006 ) ; L11 = Lovis et al. ( 2011 ) ; M09 = Mayor et al. ( 2009 ) ; M11 = Mayor et al. ( 2011 ) ; P11 = Pepe et al. ( 2011 ) ; T13a = Tuomi et al. ( 2013a ) ; T13b = Tuomi et al. ( 2013b ) .

System	Points	$\Delta t$	Model	Global rms		rms before		rms after		Remaining significant signals
HD #				std	new	std	new	std	new
10700	561	5195	0 pl.	1.23	1.12	1.13	1.06	1.89	1.43	No long-$P$ power.
	471+90									($20.0d$, $0.35m{s}^{-1}$) $\to$ ($20.0d$, $0.36m{s}^{-1}$)
										($49.5d$, $0.42m{s}^{-1}$) $\to$ ($49.6d$, $0.40m{s}^{-1}$)
										($142d$, $0.46m{s}^{-1}$) $\to$ ($142d$, $0.49m{s}^{-1}$)
										($228d$, $0.38m{s}^{-1}$) $\to$ ($230d$, $0.44m{s}^{-1}$)
										($244d$, $0.44m{s}^{-1}$) $\to$ ($243d$, $0.37m{s}^{-1}$)
										($362d$, $0.53m{s}^{-1}$) $\to$ ($360d$, $0.45m{s}^{-1}$)
										($593d$, $0.59m{s}^{-1}$) $\to$ ($596d$, $0.50m{s}^{-1}$)
										($716d$, $0.56m{s}^{-1}$) $\to$ ($723d$, $0.54m{s}^{-1}$)
										($913d$, $0.57m{s}^{-1}$) $\to$ ($948d$, $0.49m{s}^{-1}$)
										($2680d$, $0.50m{s}^{-1}$) $\to$ ($2533d$, $0.37m{s}^{-1}$)
										20.0 and 49.6 d = g and h of T13b/F17b.
										142 / 596 d = e (163 d) / f (636 d) of T13b/F17b?
										No signal at 14, 35, 94/92 d (b, c, d of T13b/F17b).
20794	594	4535	2 pl.	1.17	1.14	1.13	1.13	1.58	1.26	Small long-$P$ power.
	530+64		(P11)							(232 d, $0.53m{s}^{-1}$) $\to$ (237 d, $0.55m{s}^{-1}$);
										(617 d, $0.57m{s}^{-1}$) $\to$ (628 d, $0.49m{s}^{-1}$);
										(1877 d, $0.40m{s}^{-1}$) $\to$ (1815 d, $0.42m{s}^{-1}$).
										Do not correspond to objects of P11 or F17a.
69830	257	4817	3 pl.	1.28	1.26	1.27	1.26	1.60	1.40	Long-$P$ power.
	249+8		(L06)							(363 d, $0.68m{s}^{-1}$) $\to$ (363 d, $0.63m{s}^{-1}$)
										(4006 d, $0.58m{s}^{-1}$) $\to$ (4098 d, $0.58m{s}^{-1}$): MC
1461	266	5073	2+MC	1.68	1.64	1.65	1.60	2.54	2.84	No long-$P$ power.
	258+8		(D16)							(365, 1.27) $\to$ (366, 1.15).
154088	189	4144	1 pl.	1.78	1.74	1.78	1.74	0.70	0.03	Long-$P$ power but no significant signal.
	187+2		(M11)							The rms after, based on 2 points, is not meaningful.
			1 + d1	1.71	1.65	1.72	1.66	0.67	0.06	Very low long-$P$ power. No significant signal.
			(L11)							The rms after, based on 2 points, is not meaningful.
31527	256	5053	3 pl.	1.38	1.42	1.39	1.42	1.37	1.42	Low long-$P$ power but no signal with a
	245+11		(M11)							$p$-value¡10%.
40307	242	5059	4 + d3	1.31	1.32	1.30	1.30	1.54	1.70	No long-$P$ power. No significant signal.
	232+10		(M09)							Especially, the signal at 200 d is no more
			/T13a							significant when we fit a cubic drift
			/D16)							and no signal at 34.6 d (planet e, T13a)

Equations

$$\frac{L_{\text{gap}}}{L_{\text{nom}}}=\frac{L_{\text{est}}}{L_{\text{nom}}}-1=\frac{F_{\text{meas}}}{F_{\text{exp}}}-1$$

$$\sum _{k=0}^P{a}_k{X}_i^k={\lambda }_i$$

$$\begin{array}[]{lp{0.8\linewidth}}\{a_{k}\}_{k\in[0,P]}&{unknown} coefficients of the polynomial of the wavelength solution of order $o$. $P$ is the fixed degree of the polynomials, identical for all the orders\\ X_{i}&true position (in whatever physical unit, it be microns or nominal pixels for example) of line\leavevmode\nobreak\ $i$ on the CCD along the order ($x$-axis)\\ \lambda_{i}&(known) wavelength of line $i$ from \cite[citealias]{\@@bibref{Phrase1}{RNS14}{\@@citephrase{RNS14}}{}}\end{array}$$

$$X_i=x_i+\sum _{j=1}^G{\delta }_{i,j}{g}_j$$

$$\begin{array}[]{lp{0.8\linewidth}}x_{i}&(known) ``false'' position of line $i$ on the CCD along the order ($x$-axis) as measured in terms of pixels. Practically, it simply means that this position does not take the gaps into account.\\ G&(known) number of gaps on line $i$'s order\\ \delta_{i,j}&(known) coefficients with value 0 or 1 depending on the relative position of line $i$ with respect to gap number $j$: noting $n_{i}$ and $n_{j}$ the respective (pixel) positions of line\leavevmode\nobreak\ $i$ and gap\leavevmode\nobreak\ $j$, $\delta_{i,j}\leavevmode\nobreak\ =\leavevmode\nobreak\ \left\{\begin{array}[]{ll}1&\mathrm{if}\leavevmode\nobreak\ n_{i}>n_{j}\\ 0&\mathrm{else}\end{array}\right.$\\ g_{j}&size of gap number $j$\end{array}$$

$$\sum _{k=0}^P{a}_k{\left(x_i+\sum _{j=1}^G{\delta }_{i,j}{g}_j\right)}^k={\lambda }_i$$

$$\sum _k{a}_k{X}_i^k-\Delta {\lambda }_i={\lambda }_{i,RNS14}$$

$$\sum _k{a}_k{X}_i^k={\lambda }_{i,RNS14}$$

$$\sum _k{a}_k{X}_i^k-{\lambda }_i=0$$

$$\Delta {\lambda }_i=0$$

$$\log (\sigma \text{leavevmode}[\mathring{A}])=-0.2488\times \log (I\text{leavevmode}[{\mathrm{e}}^{-}])-1.8311$$

$$\sigma \text{leavevmode}[\mathring{A}]=0.01475\times (I\text{leavevmode}[{\mathrm{e}}^{-}]{)}^{-0.2488}$$

$$\lambda =\frac{K}{m}$$

$$\lambda =\sum _c{b}_c{\left(\frac{1}{m}\right)}^c$$

$$\sum _k{a}_k{X}^k-\sum _c{b}_c{\left(\frac{1}{m}\right)}^c=0$$

$$N\times \frac{1}{{\sigma }_{\mathrm{L}}^2}=\sum _i\frac{1}{{\sigma }_i^2}$$

$${\sigma }_{\mathrm{L}}=\sqrt{\frac{N}{\sum _i(1/{\sigma }_i^2)}}$$

Full text

*Aims. *The aim of the presented work is to improve the wavelength calibration of the HARPS spectrograph to push further its planet-detection capabilities.

*Methods. *The properties of the CCD stitching-induced pixel-size anomalies are determined with LED flat-field frames, and then a physical, gap-corrected map of the CCDs is used for the fitting model of the spectral orders. We also use a new thorium line list, based on much higher-accuracy measurements than the one used up to now. We derive new wavelength solutions for the 15 years of HARPS data, both before and after the fibre upgrade that took place in 2015.

*Results. *We demonstrate that we indeed correct the gap anomalies by computing the wavelength solutions of laser frequency comb exposures, known to have a very low dispersion, both with and without taking the gap correction into account. By comparing the rms of the most stable stars of the HARPS sample, we show that we globally decrease the radial velocity dispersion of the data, especially for the data acquired after the change of fibres of 2015. Finally, the comparative analysis of several individual systems shows that we manage to attenuate the periodogram power at one year in most cases. The analysis of the RVs derived from individual stellar lines also shows that we indeed correct the stitching-induced RV variation.

*Conclusions. *This improved calibration of the HARPS spectrograph allows to go deeper in the search for low-amplitude radial-velocity signals. This new calibration process will be further improved by combining the thorium calibration spectra with laser frequency comb and Fabry-Perot calibration spectra, and this will not only be used for HARPS but notably also for HARPS-N and the new ESPRESSO spectrograph.

1 Introduction

1.1 Context

Since Isaac Newton’s detailed study in 1666–1672 of the decomposition of light when it passes through a prism (Newton 1730), spectroscopy has been widely used for various types of studies in physics, chemistry, biology and of course astronomy. Concerning astronomy, Doppler spectroscopy, also known as the radial-velocity technique, has been used at least since 1867, when William Huggins used this technique to determine the motions of several objects with respect to the Earth (Huggins 1867, 1868). Since then, this technique has led to the discovery of many star companions and, as it was suggested by Struve in 1952, extrasolar planets. Since the first planetary-mass objects that were identified by Doppler spectroscopy, more than 20 years ago (Mayor & Queloz 1995), 786 of them have been revealed by this technique as of November 30th, 2018 (exoplanet.eu, see Schneider et al. 2011; Schneider 2015, 2016), and many others have been confirmed thanks to it. Today, the precision is already better than $1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ (Pepe et al. 2014a), what permits to already detect short-period (sub-)Earth-mass planets (see e.g. Astudillo-Defru et al. 2017), and is expected to keep improving in the coming years (see Lovis & Fischer 2010 for further historical details and Fischer et al. 2016 for the state of the art). The declared aim is to find Earth-like planets, which implies to be able to detect Earth-mass planets in the habitable zone of their stars. To reach this goal, the precision that is required is of the order of a few tens of centimetres per second or less ($9\text{\,}\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$ for the Earth around the Sun) over timescales comparable to or longer than the orbital period.

Another recent use of Doppler spectroscopy in astronomy aims at measuring several fundamental constants that rule our universe, among which the fine-structure constant and the proton-to-electron mass ratio, and determining whether they are really constant in space and time (see e.g. Murphy et al. 2003; Molaro et al. 2013a). Also, future spectrographs are expected to be able to provide a direct and model-independent determination of the expansion history of the Universe by measuring the redshift drifts of cosmologically distant sources (the Sandage-Loeb test), which requires a precision of only few $\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$ over many years (Sandage 1962; Loeb 1998; Liske et al. 2008). High-precision Doppler spectroscopy is also important for other astronomical studies such as, for example, stellar oscillations (Arentoft et al. 2008; Pepe et al. 2014b).

The High-Accuracy Radial-velocity Planet Searcher (HARPS), operational since 2003, is a high-resolution ($R=$115,000$$) fibre-fed cross-dispersed echelle spectrometer installed on the ESO 3.6-m telescope at La Silla Observatory, in Chile (Pepe et al. 2000; Mayor et al. 2003). It covers the visible spectrum from 380 to $690\text{\,}\mathrm{nm}$. As its name indicates, it is especially dedicated to the discovery and follow-up of extrasolar planets. As of 2018, it remains one of the most precise and efficient spectrographs available in the world for radial-velocity exoplanet detection, having demonstrated a long-term precision of at least $0.8\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ (Dumusque et al. 2012; Pepe et al. 2014a) and having permitted to detect radial-velocity signals with semi-amplitudes down to or even below $0.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ (Pepe et al. 2011; Dumusque et al. 2012; Feng et al. 2017a, b)111Here we use voluntarily the word “signals” instead of “planets” given that the planethood of some of these signals is uncertain or questioned, but these signals were nonetheless detected, no matter their origin.. Thorium has been used for the wavelength calibration of HARPS since the beginning, and improvements were added later. First came the Fabry-Perot interferometre, that became operational in 2011. The procedure for it, one step further than the present paper, is described in Cersullo et al. (2018). A laser-frequency comb (LFC) was added later, that had a long period of commissioning and has just been offered to the community on November 15${}^{\textrm{th}}$, 2018 (see http://www.eso.org/sci/facilities/lasilla/instruments/harps/news.html). Thanks to the last upgrades of HARPS, such as the use of a laser frequency comb (LFC) and octagonal fibres, the precision of the spectrograph is expected to become even better than $0.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ on bright stars (Lo Curto et al. 2015). However, the improvement of the instrumental performances has been a continuous process and some initial limitations still need to be overcome. Moreover, the goal of the present work is to improve all HARPS data, including the “old” ones when only the thorium lamps were available, and not only the recently acquired or future ones. These “old” data represent many years of data and should therefore not be put aside.

1.2 HARPS CCD block stitching

The HARPS CCD (full detailled description in ESO 2001) is a 4K$\times$4K CCD mosaic made of two 2K$\times$4K E2V CCDs type 44-82 with a nominal pixel size of $15\text{\,}\mathrm{\SIUnitSymbolMicro m}$. As shown on Figure 1, these individual 2K$\times$4K CCDs are themselves made of two rows of eight $1024\,\times\,512\leavevmode\nobreak\ \mathrm{pixels}$ subparts hereafter called “blocks” that arise directly from the manufacturing process of the CCD. The stitching between these blocks is not perfect, such that the pixels at the edge of each of these subparts may be slightly larger or smaller than the “central” pixels. These inter-block pixel-size anomalies are hereafter called “gaps”. The central pixels can also have small variations of size with respect to the nominal value, but these variations are expected to be much smaller than the inter-block gaps. Additionally, they are supposed to be random with zero average, and of such a low level that, with our technique, they cannot be distinguished from the pixel-to-pixel sensitivity variation222This could be achieved by, among other techniques, scanning the CCD with a laser frequency comb: see for example Wilken et al. (2010). This has been done for example by Murphy et al. (2012) and Ravi et al. (2017)., hence the fact that they will not be accounted for in our work.

The average dispersion of the HARPS spectrograph is $820\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{px}}^{-1}$ ($2.7\text{\times}{10}^{-6}$ in relative) at the middle of the orders, i.e. $0.015\text{\,}\mathrm{\SIUnitSymbolAngstrom}\text{\,}{\mathrm{px}}^{-1}$ for a wavelength of $5500\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. Earth’s yearly revolution around the Solar System barycentre (SSB) has a peak-to-peak velocity amplitude of $59.58\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, i.e. $2.0\text{\times}{10}^{-4}\text{\,}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{s}$ the speed of light in vacuum. This is equivalent to a wavelength shift of $1.1\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ at $\lambda=$5500\text{,}\mathrm{\SIUnitSymbolAngstrom}$$, i.e. to a position shift of $73\text{\,}\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}\mathrm{s}$ on HARPS’s CCD, or about one seventh of the width of a block. The 72 HARPS spectral orders are along the 512-px-side direction, so each of them crosses 7 gaps. This makes a non-negligible part of the spectral lines imaged on the CCD come and go across the stitching boundaries with a one-year periodicity, that induces a one-year signal in the final data of a number of stars observed with this instrument. This effect has been known for a long time (see e.g. Lovis & Pepe 2007, hereafter referred to as LP07, and Dumusque 2012) and was already studied by several other authors (see Wilken et al. 2010; Molaro et al. 2013b; Dumusque et al. 2015; Bauer et al. 2015).

In order to obtain an improved solar atlas, Molaro et al. (2013b) computed new wavelength solutions of HARPS spectra. They took the block-stitching effect into account by fitting separately each 512-pixel-large section of a laser frequency comb (LFC) spectrum by a third-order polynomial. They show that, by doing this, they significantly reduce the dispersion of the residuals compared to using a single third-order polynomial for the whole order assuming that all the pixels have the same size. However, one should notice that this option implies the fit of 24 coefficients for each order, that have a priori no reason to be considered as independent. This technique is only applicable for calibration sources that produce a dense enough number of lines across all CCD blocks.

Dumusque et al. (2015) analysed several stars and showed several examples for which the one-year signal is perfectly correlated with the barycentric Earth radial velocity (BERV), reinforcing the presumption that this one-year signal is related to the yearly motion of the Earth around the Sun. Moreover, a line-by-line analysis shows that, if the lines that cross the block stitching boundaries are removed from the analysis, the one-year signal disappears for some of the stars at least, while this one-year signal is further enhanced if one keeps only the gap-crossing lines. Despite the fact that the solution proposed by Dumusque et al. (2015) does not solve the problem for some of the stars, indicating that the block stitching may not be the only effect (some remnant telluric lines are for example suspected; see Cunha et al. 2014), it clearly demonstrates that the block stitching is at least part of the cause of the one-year signal.

One may argue that removing the gap-crossing lines is enough to solve the CCD-gaps part of the problem. However, it must be remarked that this technique will not solve the wavelength inaccuracy at its roots. Indeed, the position of the lines used for the wavelength solutions still assumes that all the pixels have the same size, hence keeping the slight inaccuracy in the fitted positions of these lines. Moreover, removing these lines reduces the number of usable stellar lines, thus reducing the precision that one could expect by using all the available lines, and given that these “contaminated” lines represent about 20% of all the lines (Dumusque et al. 2015), removing them has a clearly non-negligible effect on the reachable precision. The idea here is therefore to correct the thorium-line positions, in order to get a correct wavelength solution, and therefore radial velocities with an optimal precision, such that no stellar absorption line suffers from periodic systematic error and must therefore be excluded from the RV computation.

1.3 New thorium spectral line dataset

The radial-velocity technique relies on the measurement of the wavelength shift of spectral lines due to the Doppler effect and so, given that today the observed spectra are registered on CCDs, it consists fundamentally in computing the relation between a position on the CCD and the corresponding wavelength, what is called the wavelength solution. The accuracy of the latter is crucial for precise Doppler measurements. For HARPS, LP07 observed a dispersion of the residuals around the wavelength solutions reaching 50–$70\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. The origin of this high scatter is unexplained, and the knowledge of where it comes from could allow to properly correct it or at least to determine the elements that could be improved in new instruments. One can notably think about the new $10\text{\,}\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$-precision ESPRESSO spectrograph (Pepe et al. 2014b; González Hernández et al. 2017), expected to be able to detect Earth-like planets, or further planned instruments such as the HIRES instrument expected for ESO’s ELT (Maiolino et al. 2013; Zerbi et al. 2014).

Palmer & Engleman (1983), hereafter PE83, published the Los Alamos Atlas of the Thorium Spectrum, a very-high-resolution ($R=$600,000$$) catalogue of about $11\,500$ thorium spectral lines between 3000 and $11\,000\text{\,}\AA$, that had been obtained with the McMath-Pierce $1\text{\,}\mathrm{m}$ Fourier Transform Spectrometer (FTS) of the National Solar Observatory at Kitt Peak. This catalogue was used for a long time as the reference for thorium-line wavelengths and then for thorium-calibrated spectrographs (see e.g. Baranne et al. (1996) for ELODIE or LP07 for HARPS). However, even at this high resolution, the provided wavelengths are affected by systematics that lie significantly above the $10^{-10}$ relative precision aimed by sub-$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ RV spectrographs. Taking advantage of HARPS’s stability, LP07 attempted to make a new list of thorium and argon spectral lines based on HARPS ThAr spectra.

Redman et al. (2014), hereafter referred to as RNS14, made a new atlas that contains $19\,874$ lines of neutral, singly ionized and doubly ionized thorium (Th I, Th II and Th III, respectively) from 250 to $5500\text{\,}\mathrm{nm}$ ($40\,000$ to $1800\text{\,}{\mathrm{cm}}^{-1}$). They experimentally measured more than 1600 thorium lines from the near ultraviolet to the near infrared with a resolution up to $10^{7}$ with the 2 m FTS at the National Institute of Standards and Technology (NIST), and combined them with lines observed in seven other experimental studies: Giachetti et al. (1974), Zalubas & Corliss (1974), Zalubas (1976), Palmer & Engleman (1983), Engleman (2003), Lovis & Pepe (2007), and Kerber et al. (2008). RNS14 found that LP07’s thorium line list contained a large scatter with respect to their own new list, probably caused at least partially by the HARPS CCD block stitching and the non-resolution of certain blended lines.

In this work, we use therefore the new thorium line list determined by RNS14 as our wavelength reference. RNS14’s list contains two sets of wavelengths: the Ritz wavelengths, that were computed from globally optimized energy levels, and the measured wavelengths. According to RNS14, the Ritz wavelengths have a much lower uncertainty than the measured wavelengths. This is due to the fact that the Ritz wavelengths take the observed ones as input to an energy level model of the atom/ion in question, so the individual wavelengths are then effectively constrained by many observed lines. We therefore used these Ritz wavelengths as our reference.

1.4 The 2015 HARPS upgrade

In May 2015, the circular fibres of HARPS were replaced by octagonal fibres (Lo Curto et al. 2015). This change had a noticeable impact on the line profile (e.g. FWHM) and a spectral-type dependent RV offset could be observed between before and after the change of fibres. Up to now, the wavelengths used for the post-2015 data were those that had been computed by LP07 for the pre-2015 HARPS. These were therefore not optimised for this period and the wavelength solutions were consequently not as good as for the pre-2015 data. This could be seen as degraded residuals around wavelength solutions and an increase in the rms of the computed radial velocities. For these reasons, HARPS before and after the change of fibres was considered as two different instruments. The wavelength calibration improvement presented in this paper was therefore done twice independently, and the data analysis also considers the pre-2015 and post-2015 data as two separate datasets.

2 Calibration data sets

Our work uses and is based on different calibration data sets, each one for a specific task. Hereafter we describe them for a better overview.

2.1 LED flat-field frames

The first type of frames that was used is LED flat-field frames. These frames record the response of each individual pixel under a homogeneous illumination, and are then used to correct pixel-to-pixel variation effects.

Ten such frames were used for the determination of the positions and the sizes of the gaps (see Sections 3.1 and 3.2 respectively). The interest of using LED flat-field frames here is that they provide an information that is independent from the wavelength calibration frames (that use a spectral source), so a consistency check can be made after correction.

As presented in Figure 2, on these LED frames, we unsurprisingly see pixel-to-pixel brightness variations. These brightness variations are usually considered to be due to pixel-to-pixel sensitivity variations (technically referred to as Photo-Response Non-Uniformity, PRNU), assuming that all the pixels have the same size, and flat-field frames are actually usually used to correct these sensitivity variations. However, these brightness variations can also be due to effective differences of the collecting area of the pixels, for instance at the boundary between two blocks, if the electrode pattern produced by photo-lithography had not been replicated with a shift that is exactly a multiple of the pixel size. A priori, we cannot discriminate between a sensitivity difference produced by variable quantum efficiency or by a different size of the pixel. Nevertheless, a regular pattern located at the boundaries between the blocks would be a clear indication for the latter: in fact, we find every 512 pixels a column with a brightness significantly different from its neighbours. This brightness difference is much larger than the “normal” local pixel-to-pixel brightness scatter that we can see everywhere else on the frame, and the regularity of the pattern reinforces the idea that it is not due to a white noise. These columns are actually four to eight 1024-pixel-high sub-columns (four times 1024 pixels in height, and one or two pixels in width), each sub-column being darker or sometimes brighter at different levels. These are actually known since the beginning of HARPS (see e.g. ESO 2001, where a picture is shown in section 1.1.16, page 13), and have always been strongly suspected to be due (at least partially) to the gaps that we want to characterise.

2.2 Laser-frequency-comb calibration frames

A set of 82 laser-frequency-comb (LFC) calibration frames taken with HARPS between 2015-04-17T02:57:32.992 and 2015-04-17T04:20:07.352, i.e. not long before the upgrade of May 2015 (Lo Curto et al. 2015), was used to cross-check the sizes of the gaps estimated from the LED exposures. The advantage of laser frames compared to atomic-lamp frames (in our case, thorium-argon frames, see Section 2.3) is that the absolute wavelength of each line of the laser frequency comb is known with an extreme precision. Indeed, only two parameters are needed to describe the wavelengths of the whole set of laser lines: $\nu_{0}$, the anchor frequency of our laser comb, and $\Delta\nu$, the mode spacing. Then, the formula that gives the wavelength of a given line is $\nu_{n}=\nu_{0}+n\cdot\Delta\nu$, for every integer $n$. In our case, $\nu_{0}=$288.059,800\text{,}\mathrm{THz} and $\Delta\nu=$18.0\text{\,}\mathrm{GHz} (Lo Curto, private communication). Knowing the absolute wavelength of a single laser line for each order therefore gives the wavelengths of all the laser lines recorded on the CCD.

Unfortunately, the laser comb does not currently cover the whole CCD; only the orders between 26 and the last one (order 72) could then be used. This covers about three quarters of the whole CCD, especially three of the four 1024-pixel-high rows of CCD blocks, the fourth one being only covered by a single order (order 26).

2.3 Thorium-argon calibration frames

Three sets of thorium-argon calibration frames were used during this process. Those lamp frames are used as the absolute wavelength calibrators of our data. The first two sets correspond to reference calibration data of HARPS before and after its upgrade in May 2015 (Lo Curto et al. 2015), respectively. The first set is composed of 64 frames of December 2004, all with thorium-argon on both fibres, and therefore usable for the calibration of both fibres. The second one consists of 120 frames of June 2015, distributed as such: 48 frames with thorium-argon of both fibres, usable for the calibration of both fibres, and 72 with thorium-argon on fibre A and Fabry-Perot on fibre B, therefore usable only for the calibration of fibre A. Finally, the third set contains 10 frames from April 13th to April 23rd, 2015, used to compare the results with the ones obtained with the laser frames, that are of April 17th, 2015 (see Section 2.2).

The spectral lines of the first two sets were fitted and matched to our list of thorium lines, whose approximate positions on the detector are known. We got lists of lines for both fibres, both before and after the upgrade of May 2015. Before the upgrade, we counted for fibre A a list of $10\,227$ lines, and for fibre B a list of $10\,006$ lines, most of them being either thorium lines or argon lines. After the upgrade, the lists were reduced to $9944$ lines for fibre A and $9348$ lines for fibre B. The difference is due to the fact that some lines that were well fitted for before the fibre upgrade are no more well fitted after (mainly blended lines).

3 Determination of the properties of the gaps

3.1 Positions of the gaps

To determine the positions of the gaps, we used the LED flat-field frames presented in Section 2.1. A spatially regular brightness anomaly, large in intensity compared to the “normal” pixel-to-pixel brightness variation present on the rest of the image, is easily visible by eye, as shown on Figure 2. We notice that, as already stated by Molaro et al. (2013b), the positions of the large brightness anomalies are offset from exactly the multiples of 512. Molaro et al. (2013b) noted an offset of 4 pixels between the physical positions of the gaps on the CCD (at $n_{j}\leavevmode\nobreak\ =\leavevmode\nobreak\ 512\leavevmode\nobreak\ j\leavevmode\nobreak\ \mathrm{pixels}$) and the positions observed on the flat fields (at $n_{j}\leavevmode\nobreak\ =\leavevmode\nobreak\ (512\leavevmode\nobreak\ j+4)\leavevmode\nobreak\ \mathrm{pixels}$), and that 1 or 2 pixels are involved at each time. Here we find 3 or/and 4 pixels of offset, given that several of the anomalies are on two pixels of width. Indeed, the pixels that we see to be affected are 515-516, 1027-1028, 1539-1540, 2051-2052, 2563-2564, 3075-3076, and 3587-3588 (when numbering the pixels of the CCD from 1 to 4096).

3.2 Sizes of the gaps

Obviously, even if the abovementioned flat-field frames already have all individually high signal-to-noise ratios (SNR), the more frames can be accumulated, the more precise the measurement of the gap will be. We therefore stacked the ten LED flat-field frames mentioned in Section 2.1 into one “master” flat-field frame.

The first hypothesis that was made is that the gaps have a constant size between two given ($512\times 1024$)-pixel blocks. This may not be strictly true; however, the variation of the size of a gap along the border between two given blocks is expected to be very small. This assumption corresponds to assuming that the tilt between two neighbouring blocks is negligible. Practically, no significant variation in the brightness anomaly is observed along the border between two blocks, that makes us confident in the fact that this hypothesis is satisfied here.

Secondly, we assume that, by default, the pixels have a nominal size $L_{\mathrm{nom}}$, equal to $15\text{\,}\mathrm{\SIUnitSymbolMicro m}$, in both directions. The only size anomaly that we take into account is the one that we want to characterise, i.e. the gap size anomaly along the spectral orders. This assumption directly comes from the already-mentioned fact that one cannot distinguish the pixel-to-pixel intrinsic sensitivity variation from the relative pixel size anomaly, and by the fact that, as we can see on Figure 2, the amplitude of the gap-induced CCD brightness anomaly is much larger than the “normal” pixel-to-pixel brightness variation. The “size” of the pixels will therefore now refer only to their size along the direction of the spectral dispersion, and the brightness anomalies of the columns affected by the gaps are therefore considered to purely reflect their pixel-size difference along the direction of the orders, and therefore the size of the gaps.

Following the first hypothesis, we cumulated together blocks of 1024 lines of height as one single line, what increases significantly the SNR. One should notice that, assuming a random distribution of the pixel characteristics (size, sensitivity, …), this line stacking allows to significantly reduce the variation in size in the direction of the cross-dispersion, reinforcing the justification of the second hypothesis. For the same reason, this stacking also significantly reduces the “normal” pixel-to-pixel sensitivity variation in the spectral direction. At the end, we therefore got a “master flat-field frame” made of four 4096-pixel-long lines, each line containing 7 one-to-two-pixel-large “large” brightness anomalies.

The size of the gaps, relative to the nominal size of a pixel, was therefore calculated as:

[TABLE]

where $L_{\textrm{est}}$ is the estimated actual size of the affected pixel, $F_{\textrm{meas}}$ is the flux actually measured on this pixel, and $F_{\textrm{exp}}$ is the value that we would expect if the pixel had a nominal size.

One should notice that, with this convention, the negative $L_{\textrm{gap}}$ correspond to pixels smaller than nominal, i.e. to a pixel darker than “normal”, while the positive $L_{\textrm{gap}}$ correspond to pixels larger than nominal, i.e. to pixels brighter than “normal”.

The expected flux was computed by interpolating the flux of neighbouring pixels with a second-order polynomial. In each side of a given pixel $n$ for which the gap size was to be computed, the pixels selected for the interpolation were $[n-10,n-2[\cup]n+2,n+10]$.

We therefore got the estimated sizes of all the gaps for the four rows. The values that we obtained are in Table 1. For each gap, there are two values corresponding to the correction of the two pixels that are each time affected. The maximum values for individual pixels are $-34.54\text{\,}\mathrm{m}\mathrm{p}\mathrm{x}$ and +$15.27\text{\,}\mathrm{m}\mathrm{p}\mathrm{x}$, corresponding to $-28.3\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and +$12.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. Over two pixels corresponding to the same gap, the maximum values are $-62.00\text{\,}\mathrm{m}\mathrm{p}\mathrm{x}$ and +$25.76\text{\,}\mathrm{m}\mathrm{p}\mathrm{x}$, i.e. $-50.8\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and +$21.1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$.

We can compare these values with those obtained by Bauer et al. (2015). All except 4 are within 6 mpx ($\approx$5\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1}) of the values found by Bauer et al. ([2015](#bib.bib5)); most of them are within 3 mpx ($\approx$2.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}). The average and the dispersion of the differences (Bauer et al. 2015 minus our values) are respectively -0.17 mpx ($\approx$-0.14\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1}) and 6.23 mpx ($\approx$5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}).

Since we know where the spectral orders are on the CCD, we can therefore make a “pixel-size map” of the order-by-order, pixel-by-pixel correction to apply to the (extracted) spectra. This map has therefore an e2ds [extracted 2d spectrum] format, i.e. 72 lines of 4096 values. This map will be used later when we will compute the new wavelength solutions. Here is the process we used to build this map:

• •

For the orders that belong to only one of the four rows of CCD blocks (what is actually the case for all but one order), the correction is easy to find since it is simply the correction of the corresponding row.

• •

However, for the orders that cross the border between two rows, we use an averaged gap size weighted by the number of pixels extracted in each row at the given position on the order.

4 New wavelength solutions

In an analogous way as LP07, the idea was to fit the wavelengths of the spectral lines as a polynomial of their positions on the CCD. The difference between LP07 and the present work is that here the gap sizes are taken into account in the equations and that we use the new thorium wavelengths of RNS14.

4.1 Mathematical formula

The basic idea is to keep LP07’s general formula. For a given order $o$, and a given line $i$ of this order, the equation is:

[TABLE]

where:

[TABLE]

The $X_{i}$ of this formula is supposed to be the true position of the line on the CCD, actually corresponding to the nominal-pixel-scaled pixel-number position plus the gaps:

[TABLE]

where:

[TABLE]

The $\{n_{j}\}_{j\in[1,G]}$ and the $\{g_{j}\}_{j\in[1,G]}$ are those determined respectively in Sections 3.1 and  3.2.

When we insert Formula (3) inside Formula (2), we get:

[TABLE]

The only unknowns are the $\{a_{k}\}_{k\in[0,P]}$, and equation (4) is then linear with respect to the unknowns. For simplicity, we will from now on use the $X_{i}$ (= $x_{i}+\sum\limits_{j=1}^{G}\delta_{i,j}g_{j})$ in the writing of the equations.

Practically, equation (4) will be different depending on the type of line that we are dealing with. We have 4 types of lines: the thorium lines that appear on two different orders (“duplicate thorium lines”), the thorium lines that appear on only one order (“unique thorium lines”), the non-thorium lines that appear on two orders (“duplicate non-thorium lines”) and the non-thorium lines that appear on only one order (“unique non-thorium lines”).

Duplicate thorium lines:

[TABLE] 2. 2.

Unique thorium lines:

[TABLE] 3. 3.

Duplicate non-thorium lines:

[TABLE] 4. 4.

Unique non-thorium lines: no equation

In these equations, the $\Delta\lambda_{i}$ is the unknown systematic offset of line $i$ (for thorium lines), $\lambda_{i,\;{\rm RNS14}}$ is the RNS14 wavelength of line $i$ (for thorium lines) and $\lambda_{i}$ is the unknown wavelength of line $i$ (for unidentified or argon lines). Argon lines are included in the group of lines with unknown wavelengths because they do not appear in RNS14’s list and because they are significantly more sensitive to physical conditions in the lamp than the thorium lines (see LP07). As in LP07, the $\Delta\lambda_{i}$ is the same for both occurrences of a given duplicate thorium line, since systematic effects should affect in the same way these two occurrences. Coupling the orders in this ways also strengthen the reliability of the wavelength solutions towards the order edges, where the wavelength solutions are usually more poorly constrained.

4.2 Thorium-line systematics

The individual systematics $\Delta\lambda_{i}$ are unknowns to be determined. To understand their origin, we computed the wavelength solutions for each individual order based on RNS14’s wavelengths. Only thorium lines are therefore used, both unique and duplicate, given that only thorium lines are present in RNS14’s catalogue. Since the orders are decoupled, the system contains only equations of type (6). The residuals of the two occurrences of the duplicate thorium lines (the left-hand-side occurrence being the one of the order following the one of the right-hand-side occurrence) are plotted against one another in Fig. 3, with a superimposed 1-to-1 line. We clearly see that there is a good 1-to-1 correlation between the left-hand-side residuals and the right-hand-side ones. Moreover, we can see that the distribution of the residuals is centered on zero. From this observation, two things can be deduced: first, the residuals (so basically the $\Delta\lambda_{i}$) are on average equal to zero; second, the $\Delta\lambda_{i}$ of the equations corresponding to a same duplicate line should all be considered to be equal. In addition, it means that the $\Delta\lambda_{i}$ represent systematics that are not related to the CCD, the optics or the PSF, but rather are intrinsic to the thorium spectrum; they could for example be related to blends that are not resolved by HARPS. These elements will add significant constraints on the $\Delta\lambda_{i}$. These global properties of the $\Delta\lambda_{i}$ translate into these additional equations (one for each $\Delta\lambda_{i}$, i.e. one for each duplicate thorium line):

[TABLE]

This has to be understood in a least-squares sense: these equations will ensure that, in the global solution, the properly normalised (see below) $\Delta\lambda_{i}$ will have a mean value close to zero and a variance close to one. The weight of these equations should be representative of the dispersion of the corresponding residuals (cf. Section 4.3).

4.3 Weighting of the equations

We need to assign the proper weight to each equation in order to properly solve our weighted linear least-squares problem. Uncertainties for each equation will be computed differently depending on the equation type. For the equations of types (5) and (7), uncertainties are the quadratic sum of photon noise (considered here as the wavelength-converted uncertainty on the fitted position of the line) and a fixed detector-related noise (0.01 pixel). This detector-related noise is a noise floor such that we cannot measure the position of a line with an accuracy better than this value. This value of 0.01 px was chosen because, as we will see in Section 5.1.2 and Figure 5, the residuals of the laser wavelength solutions have a rms of about this value, whose origin is most likely instrumental (e.g. intra-pixel sensitivity variations) and not from the (much lower) photon noise. For the equations of type (6), uncertainties are given by the quadratic sum of the uncertainty on the systematic offset $\Delta\lambda_{i}$, the line-position uncertainty and RNS14’s wavelength uncertainty. Finally, equations of type (8) are weighted by the quadratic sum of the uncertainty on the systematic offset and the RNS14 wavelength uncertainty.

To determine the true uncertainty of each line (for single lines) or couple of lines (for actually duplicate lines), we can go back to Fig. 3. In addition to comparing the left-hand-side residuals to the right-hand-side ones of the duplicate thorium lines, this figure compares these residuals to the actual intensity (flux) of these lines. One has to note that, here, the flux are the flux measured by HARPS, corrected from the blaze and average over the two occurrences. One can see that the dispersion of the residuals seems to be lower for lines with higher intensities. More specifically, if we plot the rms of the residuals with respect to the flux, we find a linear relation between the logarithm of both quantities. This is shown on Fig. 4, where we did 10 bins of equal size on the log(flux) scale. This correlation allows us to use this flux-based fitted rms as the uncertainty on the systematic offset $\Delta\lambda_{i}$ mentioned before. This makes sense given that only thorium lines have equations using this uncertainty (equations of types (6) and (8)).

We fitted a linear relation between log(rms) and log(flux) and we therefore get the following relation between the uncertainty on the systematic offset of the thorium lines and their flux:

[TABLE]

i.e.

[TABLE]

For a line at $5500\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, this relation corresponds to an uncertainty on the systematic offset equivalent to $25\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for a flux of 1 million electrons and to about $150\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for a line with only 1000 electrons. This is similar to what LP07 had found for their flux-dispersion relationship, where they got uncertainty between 20 and $160\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. This shows that the possible errors in the wavelengths of PE83, assumed by LP07 to be the cause of these systematic corrections, can only be part of it and that necessarily additional causes are also responsible of these systematics. One can think of line blends due to the lower resolution of HARPS compared to RNS14 or other systematic effects.

4.4 Littrow test

From the spectrograph optical design, the wavelength at a given position along the main dispersion must follow a smooth relation from an order to the next. From the grating equation, we can write that, for a given $x$, this relation is:

[TABLE]

where $K$ is a constant and $m$ is the physical number of the echelle order (161-89 except 115 for fibre A; 161-89 except 116 and 115 for fibre B). In order to include optical aberrations, we decided to actually use a third-degree polynomial:

[TABLE]

For a given $x$ and a given order, we know that this point must satisfy both the wavelength solution of its order (equation 2) and the Littrow equation of its position (equation 12):

[TABLE]

In order to constrain more strongly the wavelength solutions, we add to the system of equations these Littrow equations for all the orders at five different positions (the same five positions for all the orders). We gave the same weight $\frac{1}{\sigma_{\mathrm{L}}^{2}}$ to these $N$ equations. The cumulated weight of these equations was set to be equal to the total weight coming from the addition of the weights corresponding to the uncertainty on the position of all the lines $i$ used:

[TABLE]

The corresponding uncertainty attributed to these equations is therefore:

[TABLE]

We chose to give the same weight to each of these equations because they are not directly linked to the line content of the spectrum.There is therefore a priori no reason to favor one position over another. Also, the overall weight of these equations was set such that it does not “hide” the other equations, but also such that it has a noticeable influence. From a “mechanical” point of view, there is no reason that one direction (along the ordres) be privileged with respect to the other (across the ordres). Equating the overall weight of these two types of equations therefore reflects this.

4.5 New wavelength solutions

Before solving the linear least-squares problem, we first built a master line list from the set of thorium frames described in Section 2.3. Measured line positions are averaged over a large number of frames to minimise photon-noise errors on line positions. The system of equations was then solved a first time to identify outliers, and then a second iteration was done with only the lines that were less than $7\leavevmode\nobreak\ \sigma$ away from the solution of the first run. This is the solution of this second run that was used for the final attribution of wavelengths. For the dataset of December 2004, we get a final reduced chi-square of 1.02 for fibre A and 0.81 for fibre B. For the dataset of June 2015, we get a final chi-square of 1.42 for fibre A and 1.07 for fibre B.

The new wavelengths obtained from this procedure are then saved in a new reference line list that is used to reprocess all individual wavelength calibration exposures acquired with HARPS over 15 years. This new list is specific to HARPS and contains the updated wavelengths that include the line systematics $\Delta\lambda_{i}$. The typical dispersion of the residuals around the new wavelength solutions is $10\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for the pre-upgrade data and $15\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for the post-upgrade data. While the pre-upgrade solutions show similar dispersions compared to the old reduction, we observe a significant improvement for the post-upgrade solutions (from 30 to $15\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$).

5 Results

5.1 Test of the gap correction with the laser frequency comb

5.1.1 Principle

We used laser-frequency-comb frames to check our gap size estimates. The advantage of laser frames is that the absolute wavelength of each line is known with an extreme precision. When we do the wavelength solution by fitting a polynomial $P$ such that $\lambda=P(x)$, where $\lambda$ is the wavelength and $x$ is the position in number of pixels, the gaps should cause jumps in the residuals, given that “pixel number” is not a physical unit of length (given the pixel-size variations). On the contrary, if we take into account the gaps by using the actual physical position of the lines as the $x$ coordinate, no jump should remain. Obviously, the gaps would be visible only if their effect is larger than the dispersion of the residuals; however, for the laser, the dispersion is expected to be very small, hopefully smaller than the effect of the gaps (the value of this dispersion is checked afterwards). This is why we used these laser frames to check the estimated values and not the thorium frames, that are known to have a much larger dispersion, about 50–$70\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, larger than the RV-converted size of the gaps. The set of 82 laser frames mentioned above (Section 2.2) was therefore used to perform this test.

As mentioned above, the laser comb does not currently cover the whole CCD, but only the orders from 25 to 69. Three of the four 1024-pixel-high rows of CCD blocks are fully covered, while only one single order is on the fourth one. This could limit the consistency of our verification; however, it seems reasonable to assume that, if the estimated sizes of the gaps look consistent for all the orders of the three fully-covered rows, and that the single order of the fourth raw also gives consistent results, the full set of size estimates is consistent with the actual sizes of the gaps. This is therefore what we considered in the rest of this study.

5.1.2 Results

The position of the lines were first fitted in pixels, i.e. without correcting the gaps, and then the wavelength solutions were computed. Clear jumps were visible in the residuals. The same was then done using a corrected position of the lines, corresponding to the pixel-position corrected by the gaps we estimated before. If these corrected positions correspond to the actual physical position of the lines, we should not see any jump in the residuals. As can be seen on the spectral orders shown in Figure 5, by doing this correction, we do not see any noticeable jump, what indicates that the estimated sizes of the gaps are correct within the residual dispersion of about $7\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$.

This low residual dispersion of the laser solutions, about $7\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, i.e. about $9\text{\,}\mathrm{mpx}$, indicates that the intrinsic dispersion due to instrumental effects is not larger than this value. This rules out this possibility as a major explanation of the larger dispersion (~50–$70\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$) observed with thorium frames. It reinforces the hints that the large scatter seen on thorium solutions is mainly due to intrinsic features of the thorium spectrum.

5.2 Comparison of the thorium and laser solutions

We further compare the laser and thorium wavelength solutions. In Figure 6, we plot the difference between the wavelength solutions from the laser and the thorium, for both the uncorrected (standard) and corrected (new) thorium solutions. As can be seen, the difference in much flatter in the case of the new version of the data reduction software (hereafter DRS) than with the standard one. In the case of the standard reduction, we can also clearly see the jumps where the gaps are, while these are not present with the new reduction.

One can however notice that the new wavelength solution has an average offset of about $43\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ with respect to the laser solution, while the standard thorium solution has an offset of about $49\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. The offset seems to be larger for longer wavelengths. As we explained in Section 2.2, the theoretical wavelengths of the laser rfequency comb (LFC) respect the formula $\nu_{n}=\nu_{0}+n\cdot\Delta\nu$, where $\Delta\nu$ ($18\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ in our case) is known at the Hz level and $\nu_{0}$ was set at $288.0598\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ (LFC team, private communication). We note that a fixed error $\Delta\nu_{0}$ on $\nu_{0}$ would translate into an error in RV proportional to the wavelength: $\Delta v=\Delta\nu_{0}\cdot\lambda$. An error $\Delta\nu_{0}=$100\text{,}\mathrm{M}\mathrm{H}\mathrm{z}$$ would translate into an offset $\Delta v=40$, 50, and $60\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for a wavelength $\lambda=4000$, 5000, and $6000\text{\,}\AA$ respectively. This is close to what we observe. The LFC frames we are using were obtained during the LFC commissioning, during which many changes to the comb properties occurred. In particular, $\nu_{0}$ was varied by a few hundred MHz over the course of the commissioning. A mistake in the reported $\nu_{0}$ for the specific frames we are using is thus a plausible explanation.

This is further strengthened by the fact that Wilken et al. (2010) did not observe any offset between LFC and thorium on the single spectral order that was available at that time. While the thorium solutions have not changed in time, the anchor frequency $\nu_{0}$ of the LFC has over the different commissioning periods. Unfortunately, we are forced to use the 2015 LFC data because these are the ones providing by far the largest spectral coverage.

We believe that the other possible causes for the offset are less likely than a $\nu_{0}$ problem. Thorium wavelengths from RNS14 have been remeasured several times by the NIST teams over the years, yielding consistent results. Furthermore, any systematic effects in fitting LFC and thorium line positions should also have appeared in Wilken et al. (2010), which is not the case. In conclusion, a wrong $\nu_{0}$ is the most likely culprit for the RV offset we observe. This is unfortunate but we have no means to investigate this further. We stress that, if our hypothesis is correct, the thorium solutions will be correct also in absolute terms.

We can note that, if we compare the two thorium wavelength solutions, we have a mean offset of $c\times\frac{\lambda_{\mathrm{new}}-\lambda_{\mathrm{old}}}{\lambda_{\mathrm{laser}}}\approx$-6\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1}. This is fully compatible with the correction factor $\kappa=0.22\times 10^{-7}$ applied by [RNS14](#bib.bib51) on [PE83](#bib.bib45)’s wavenumbers, that corresponds to an offset of $c\times\frac{\lambda_{\mathrm{RNS14}}-\lambda_{\mathrm{PE83}}}{\lambda_{\mathrm{RNS14}}}=$-6.6\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}. This has to be expected since the recomputation of the wavelengths by LP07 and by us are fit, for the absolute values, on PE83 and RNS14 respectively.

5.3 Global results on stars

A global test to check whether our implementation improved the velocities is to check whether the rms of the most stable stars has decreased with the new DRS with respect to the standard DRS. This comparison is done on Fig. 7 for the observations before the change of the fibres (observations until May 2015 included, Fig. 7a) and after the change of fibres (since June 2015, Fig. 7b). In both cases, we took the subset of objects observed at least 10 times during the period by HARPS that have a radial velocity rms below $2\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ with both the standard and the new versions of the DRS. We therefore have a set of 124 stars for before the change of fibres and 21 stars for after the change of fibres.

For the observations before May 2015, we can see from the cumulative histograms that the new DRS provides a clear improvement for the most stable stars. Indeed, the 10 stars with the lowest RV dispersions with both reductions all have their rms improved with the new DRS. For the observations after the change of the fibres, we can see that the rms are globally significantly improved. The 20 most stable stars all have their rms improved (except one), in many cases by a significant margin. An example of this improvement is shown on Figure 8, where we can see the RVs of star HD 220507 after the change of fibre for both the standard and the new versions of the DRS. For this star, where no significant signal is observed, the rms on 82 points covering 803 days after the change of fibres evolves from $1.70\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ to $1.23\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. This improvement is especially due to the recomputation of the thorium wavelengths for the post-upgrade HARPS as explained in Section 4.

5.4 Analysis of individual systems

Given the consistency of the previous results, we further test the new reduction on a set of well-observed stars followed over many years with HARPS.

All the data obtained with HARPS for these stars, spanning up to nearly 5200 days (more than $14\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$) in total, were reduced with the new pipeline taking into account the block stitching and the newly-defined thorium wavelengths. These results were compared to the ones obtained with the standard pipeline, that did not take into account the gaps and used the wavelengths as defined in LP07. Notably, we compared the periodograms, to check that the results concerning the known planetary signals and the long-term trend remained basically unchanged, and that at the same time the previously-present 1-year-period signal and its aliases had disappeared, or at least significantly decreased. The periodograms and Keplerian fits presented below were obtained thanks to the Data & Analysis Center for Exoplanets333https://dace.unige.ch (DACE; see Buchschacher et al. 2015). Keplerian model initial conditions are computed using the formalism described in Delisle et al. (2016).

The systems that we analysed are HD10700, HD20794, HD69830, HD1461, HD154088, HD31527, and HD~40307.

For all the systems, a noise level of $0.75\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ was added quadratically to the DRS-computed RV uncertainties to take into account stellar jitter and the typical instrumental noise floor seen in HARPS data. In addition to the signals fitted, we allow an offset to be adjusted between the two sets of data before and after the change of fibres. All the false-alarm probabilities (hereafter FAP) given here were computed analytically by DACE following Baluev (2008).

In Table 2 are listed the planetary systems that we studied, the number of points that we used, the timespan covered by the data, the model that we fitted, the rms of the data after removing the aforementioned model (global rms, rms before and rms after the change of fibres, with both the standard and the new versions of the DRS), as well as complementary comments on the remaining signals. The detailed model for each system is explained on individuals sections.

5.4.1 HD 10700

HD 10700 (Tau Ceti) is a G8V star located 3.7 parsecs away from the Solar System. It is one of the most stable stars in radial velocities with a long-term dispersion of \sim$$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1} (Pepe et al. 2011; Tuomi et al. 2013b). Both with the standard and the new versions of the DRS, we find a lot of peaks in the periodogram with a FAP below 1%. The 1-year signal becomes slightly less significant with the new version (FAP of 0.49%) compared to the standard one (FAP of 0.057%). While the peaks around 600 and 900 days have a very high significance with the standard version, they have a lower one with the new version, while still having a FAP well below 0.1%. On the contrary, the peak around 140 days becomes significantly stronger. The 20-day signal, that corresponds to planet g of Tuomi et al. (2013b) and Feng et al. (2017b), is basically undetected in the periodogram with the standard version while it becomes marginally significant with the new one. The signal at 49.6 days, that corresponds to planet h, is detected in both cases with a FAP below 1%. There seems to be a signal around 230-240 days, even if its double-peak structure makes it unclear, especially given that one of the two peaks is stronger with one version while the other peak is stronger with the other version. A marginally-significant peak is detected around 720 days with the standard version while it is not significant at all with the new version. Finally, the long-period signal around 2700 days, that was very significant with the standard version, becomes insignificant with the new one. While two of these signals (20.0 and 49.6 days) seem to correspond to signals detected by Tuomi et al. (2013b) and Feng et al. (2017b), the others do not seem to correspond. We may tentatively associate the signal around 600 days to their planet f and the signal around 142 days to planet e, but the periods are quite different and this correspondence looks doubtful. Concerning the other planets (b, c and d), no signal is detected at the corresponding announced periods (respectively 14, 35 and 94/92 days).

As a conclusion, we urge caution in the interpretation of the many weak RV signals seen in HD 10700, as we demonstrate that their amplitude and significance depend on the pipeline version used and the data analysis method.

5.4.2 HD 20794

HD 20794 is a G6V star located 6.0 parsecs away from the Solar System. Two signals, corresponding to planets HD 20794 b and HD 20794 d respectively (Pepe et al. 2011), are detected in our periodogram analysis, at $18.317\pm 0.005\text{\,}\mathrm{d}$ and $89.65\pm 0.10\text{\,}\mathrm{d}$ with the standard DRS and at $18.320\pm 0.005\text{\,}\mathrm{d}$ and $89.76\pm 0.12\text{\,}\mathrm{d}$ with the new DRS. For the amplitudes of these signals, we find respectively $(0.67\pm 0.07)$ and $0.71\pm 0.07\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for HD 20794 b, and $(0.74\pm 0.09)$ and $0.57\pm 0.07\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for HD 20794 d. While the difference is within 1 $\sigma$ for planet b, it is more noticeable for planet d, 22% and nearly 2 $\sigma$ lower with the new DRS than with the standard one. The fitted eccentricities are also quite lower with the new DRS, changing from $0.18\pm 0.10$ (b) and $0.45\pm 0.09$ (d) with the standard version of the DRS to $0.11\pm 0.10$ (b) and $0.29\pm 0.11$ (d). The other candidate planets, HD 20794 g, c, e and f at respectively 12, 43, 147 and 330 days (Pepe et al. 2011; Feng et al. 2017a), are not found. However, we see signals around 235, 620 and 1900 days. The first two of these signals may correspond to those around 250 and 680 days found by Feng et al. (2017a), that they interpret respectively as a yearly alias of their signal at 147 days and a harmonic of their signal at 330 days. Given that we do not detect the signals at 147 and 330 days, the actual situation may be the opposite. Again, caution is in order when interpreting the nature of these additional signals, whose properties vary according to the pipeline version used. A full study is beyond the scope of this paper.

5.4.3 HD 69830

HD 69830 is a K0V star located 12.6 parsecs away from the Solar System. The three planets announced in Lovis et al. (2006), with periods of 8.67, 31.6 and about 200 days, are recovered. After fitting these three planets, two peaks remain visible below the significance threshold: one around one year and a long-period one around 4000 days, that may correspond to the magnetic cycle of the star. The one-year peak is slightly less significant with the new DRS while the long-period peak is slightly more significant. Since the long-period signal is not significant, we did not fit it here. However, we checked and found that fitting it ($P\approx$3680\text{,}\mathrm{d}$,K\approx$1.1\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1}$,e\approx 0.8$) does not change much the values we obtain for the planets, all remaining well within one sigma. For the innermost planet, we obtain a period of $8.6682\pm 0.0004\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $3.47\pm 0.11\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.07\pm 0.03\text{\,}$. For the second planet, the period is $31.601\pm 0.006\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, the amplitude $2.66\pm 0.13\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and the eccentricity $0.06\pm 0.05\text{\,}$. Finally, the outermost planet is found to have a period of $200.32\pm 0.36\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $1.85\pm 0.12\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.12\pm 0.06\text{\,}$ with the standard version of the DRS while we find a period of $200.57\pm 0.37\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $1.77\pm 0.12\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.10\pm 0.06\text{\,}$ with the new version.

5.4.4 HD 1461

HD 1461 is a G0V star located 23.4 parsecs away from the Solar System. Díaz et al. (2016) announced two planets (periods of 5.77 and $13.5\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$) around this star as well as a magnetic cycle (period of $3500\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$). We do find these three signals with both the standard and the new version of the DRS. Without fitting the magnetic cycle, we find a signal at 37.3 days with both methods, but while its FAP is slightly below 1% with the standard version of the DRS, it is slightly above with the new version. This signal disappears once we fit the magnetic cycle. Both with and without fitting the magnetic cycle, we find a significant signal around 1 year, even if we can notice that the amplitude is a bit smaller with the new version of the DRS.

All the parameter values obtained for this system from the old and new versions of the DRS are well within $1\leavevmode\nobreak\ \sigma$ of each other, so here we give only the values from the new DRS. The fit of the magnetic cycle gives a period of $3582\pm 80\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $3.28\pm 0.28\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.51\pm 0.07\text{\,}$. Concerning the innermost planet, we find a period of $5.7715\pm 0.0004\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $2.26\pm 0.15\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.11\pm 0.07\text{\,}$. For the second planet, we get a period of $13.506\pm 0.003\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $1.60\pm 0.16\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.15\pm 0.10\text{\,}$.

5.4.5 HD 154088

HD 154088 is a K0IV star located 17.8 parsecs away from the Solar System. The single planet announced by Mayor et al. (2011) at a period of 18.6 days around this star is easily recovered. No other significant signal is found in the HARPS RV data of this star, even if we can notice the presence of some low power at long periods in the periodogram. We therefore decided to model these data both with a single Keplerian and with a Keplerian plus a linear drift. In both cases, the values from the old and the new DRS are all compatible within $1\leavevmode\nobreak\ \sigma$ so, as before, we give only the values from the new DRS.

From the model with only a Keplerian, we find a period of $18.556\pm 0.009\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, nearly $2\leavevmode\nobreak\ \sigma$ lower than the $18.596\pm 0.021\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ in Mayor et al. (2011). The amplitude is found to be $1.68\pm 0.19\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, slightly lower than but well within $1\leavevmode\nobreak\ \sigma$ of the $1.78\pm 0.31\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ of Mayor et al. (2011). We find an eccentricity of $0.10\pm 0.11\text{\,}$, noticeably lower than the $0.38\pm 0.15\text{\,}$ in Mayor et al. (2011).

According to Lovis et al. (2011), a long-term magnetic cycle is present in this star, but its period is very unconstrained, given as $8000^{+\infty}_{-6487}$ days in this paper. Here, we can clearly notice a trend in the data. Fitting a linear drift significantly reduces the dispersion of the residuals. This drift has an amplitude of $0.33\pm 0.07\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{yr}}^{-1}$. In this case, we find for the planet a period of $18.568\pm 0.008\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, slightly larger than before but still lower than what had found Mayor et al. (2011). The amplitude is found to be $1.64\pm 0.23\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, slightly lower than before but well within $1\leavevmode\nobreak\ \sigma$. We find an eccentricity of $0.28\pm 0.12\text{\,}$, noticeably larger than before but still lower than what had found Mayor et al. (2011). Even if its significance is extremely low in all cases, one can notice that the peak around 360 days is smaller with the new DRS than with the standard one.

5.4.6 HD 31527

HD 31527 is a G0V star located 38.4 parsecs away from the Solar System. Mayor et al. (2011) announced the detection of three planets around this star, at periods of 16, 51 and 274 days. All three signals are easily recovered with both versions of the DRS and no other significant signal is found in the HARPS RV data of this star. The highest remaining signal is a multi-peak between 300-400 days, but with negligible significance.

Both versions of the DRS give values compatible within $1\leavevmode\nobreak\ \sigma$, except for the amplitudes of planet d, that are nonetheless within less than $2\leavevmode\nobreak\ \sigma$. Therefore, here again, we give only the values from the new DRS. Concerning the planets, HD 31527 b is found to have a period of $16.557\pm 0.003\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $2.84\pm 0.13\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.16\pm 0.05\text{\,}$. HD 31527 c has a period of $51.20\pm 0.03\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $2.70\pm 0.14\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.08\pm 0.05\text{\,}$. Finally, HD 31527 d has a period of $274.2\pm 1.0\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $1.77\pm 0.36\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.67\pm 0.08\text{\,}$. All these values are compatible with those given by Mayor et al. (2011), even if we find a noticeably larger value for the eccentricity of planet d. One can also note that, with 89 more points and nearly twice as long a timespan as in Mayor et al. (2011), we have significantly lower uncertainties on the parameters of the outermost planet.

5.4.7 HD 40307

HD 40307 is a K3V star located 13 parsecs away from the Solar System. Mayor et al. (2009) discovered three planets (b, c and d) around this star, at periods of 4.3, 9.6 and 20.5 days, that we easily recover. Tuomi et al. (2013a) announced the detection of three more planets (e, f and g), at periods of 34.6, 51.8 and 197.8 days. Díaz et al. (2016) confirmed the three planets of Mayor et al. (2009) as well as the 51.6-day planet of Tuomi et al. (2013a). While they find a signal at 205 days that may correspond to Tuomi et al. (2013a)’s planet g, they conclude that its significance and interpretation are doubtful. They also find that the six-planet model including the planet at 34.6 days is strongly disfavoured. Here, we recover the three planets of Mayor et al. (2009) as well as the signal at 51.6 days. As found by Díaz et al. (2016), we detect a long-term signal, probably the magnetic cycle of the star, that we fitted by a third-degree polynomial. Once this is done, no significant signal remain. Especially, the signal around 200 days has a very low significance and no signal is found around 34.6 days.

The two DRS give values compatible within $1\leavevmode\nobreak\ \sigma$, so here we also give only the values from the new version. The innermost planet has a period of $4.3114\pm 0.0002\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $1.84\pm 0.14\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.12\pm 0.07\text{\,}$. The second planet has a period of $9.6210\pm 0.0008\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $2.29\pm 0.13\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.05\pm 0.06\text{\,}$. The third planet has a period of $20.412\pm 0.004\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $2.31\pm 0.14\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.07\pm 0.05\text{\,}$. The fourth signal has a period of $51.68\pm 0.07\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$, an amplitude of $0.78\pm 0.13\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and an eccentricity of $0.19\pm 0.16\text{\,}$. Finally, the drift parameters are $0.60\pm 0.07\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{yr}}^{-1}$, $-0.13\pm 0.02\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{yr}}^{-2}$ and $-0.02\pm 0.00\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{yr}}^{-3}$.

5.4.8 Conclusion

We can see that the global rms is generally slightly smaller with the new version than with the standard version. For the subset of data obtained after the change of fibres, the rms is generally significantly better, as shown by the results of Section 5.3 and in Figure 8.

Concerning the fitted signals (confirmed planets and magnetic cycles), all the amplitudes are compatible within 2 $\sigma$, most of them within 1 $\sigma$, between the two versions of the DRS, as well as between these results and those in the literature.

We also compare the periodogram power around 1 year between the two versions of the DRS. Among the 7 systems studied, the periodogram power at 1 year decreases with the new reduction in 6 cases. In the only case where it slightly increases (HD 31527), the power at 1 year remains below the significance threshold. This demonstrates that spurious instrument-related 1-year signals are significantly reduced with the new DRS. Overall, we note that several systems exhibit a complex periodogram, with significant signals not previously identified as planets at various periods. Through the window function of the data, this excess power will partially leak into other frequencies, in particular at one year, which further complicates our analysis. Nevertheless, among the 7 systems, only one (HD 1461) seems to have a significant signal at 1 year with the new reduction. It is possible that this particular signal is in fact of astrophysical origin.

5.5 Test on individual stellar spectral lines

To see the correction that our new wavelength solution brings at the location of the stitchings boundaries, we derived the velocity of a few spectral lines of the star HD 1461 passing over CCD stitching boundaries using the method described in Dumusque (2018). In Fig. 17, we show the RV of a few spectral lines passing over CCD stitchings boundaries derived using the standard data reduction software (DRS), in which the wavelength solution is estimated using a unique pixel size, and derived using the formalism described in this paper, in which the wavelength solution is obtained using the real pixel size. We also show in the bottom panels of Fig. 17 the difference between those two sets of RV.

As we can see in the standard DRS RVs, there is a significant drift when a spectral line goes over a CCD stitching boundary, shown in Fig. 17 by red vertical dashed lines. In addition, by looking at the blue horizontal lines that represent the median RV before and after the spectral line passes over the CCD stitching boundaries, we see that these two median RVs are different. This can be explained by the fact that because the pixel at the stitching does not have the same size as the other pixels, the wavelength solution after the line passed over the stitching boundary is under or over-estimated, thus inducing an offset in RV. When looking at the RVs derived using the method described in this paper, the offset between the median RV before and after the stitching boundary disappears, proving that our new wavelength solution considers correctly the different pixel sizes at the stitching positions. This reinforces the idea that the remaining signal around one year for this specific star is not due to this stitching effect but could well be of astrophysical origin.

6 Conclusion

We manage to estimate the size of the gaps of the CCDs of HARPS that induce a one-year-period parasitic signal observed in the periodogram of several stars. These sizes range from $-0.034$ to $+0.015\text{\,}\mathrm{px}$. The validity of this correction was verified with laser frequency comb exposures. 2. 2.

We reprocess all the HARPS wavelength calibration exposures with an improved model taking into account the gaps and the most recent thorium-line wavelengths. We verified that the new wavelength solutions show smaller wavelength distortions than before by comparing to laser-frequency-comb data. Moreover, the new solution are smooth accross bolck stitching boundaries, while the old solutions show significant jumps. 3. 3.

While the RV dispersion of the data obtained before the change of fibres of May 2015 is slightly but consistantly improved, the rms of the data obtained after the change of fibres is significantly better. 4. 4.

A reanalysis of a number of planetary systems with the new pipeline is able to recover all the clearly significant planetary and magnetic signals published in the literature, with fitted parameters fully compatible with the published ones. On the contrary, several weak signals close to the noise level have different properties depending on the pipeline version used. 5. 5.

When we use the new DRS, we systematically reduce the periodogram power at a period of one year in the analysed systems. 6. 6.

The analysis of the RVs derived from individual stellar lines shows that we indeed correct the RV systematic offsets induced by the crossing of stitching boundaries.

The same wavelength recalibration process can be applied to any other current and future spectrographs, such as the very-high-precision instruments HARPS-N, ESPRESSO, NIRPS, SPIRou and ELT-HIRES.

The next step to further improve the wavelength calibration of fibre-fed spectrographs will consist in the combination of Fabry-Perot and thorium exposures (see e.g. Wildi et al. 2011; Schäfer & Reiners 2012; Halverson et al. 2014; Schwab et al. 2015; Bauer et al. 2015; Cersullo et al. 2017), coupled with the improved procedure described here. This is beyond the scope of the present paper and will be described in details in Cersullo et al. (2018).

Acknowledgements.

We thank the anonymous referee for his/her remarks that helped a lot improve the article. We acknowledge the support of the Swiss National Science Foundation. X. D. thanks the Branco Weiss Fellowship–Society in Science for its financial support. This publication makes use of the Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. The NCCR are a research instrument of the Swiss National Science Foundation.