© ESO 2019

1 Introduction

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 (1952), extrasolar planets. Since the first planetary-mass objects that were identified by Doppler spectroscopy more than 20 yr ago (Mayor & Queloz 1995), 786 such objects have been revealed by this technique as of November 30th, 2018¹ (see Schneider et al. 2011; Schneider 2015, 2016), and many others have been confirmed thanks to it. Today, the precision is already greater than 1 m s⁻¹ (Pepe et al. 2014a), already permitting us to detect short-period (sub-)Earth-mass planets (see e.g. Astudillo-Defru et al. 2017). This precision 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). One declared aim of exoplanet research is to find Earth-like planets, which implies that we are 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 cm s⁻¹ 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 to measure several fundamental constants that rule our universe, including 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 cm s⁻¹ over many years (Sandage 1962; Loeb 1998; Liske et al. 2008). High-precision Doppler spectroscopy is also important for other astronomical studies such as, 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 nm. As its name indicates, it is especially dedicated tothe 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 m s⁻¹ (Dumusque et al. 2012; Pepe et al. 2014a) and having permitted the detection of radial-velocity signals with semi-amplitudes down to or even below 0.5 m s⁻¹ (Pepe et al. 2011; Dumusque et al. 2012; Feng et al. 2017a,b)². Thorium has been used for the wavelength calibration of HARPS since the beginning, and improvements were added later.First came the Fabry–Perot interferometre, which became operational in 2011. The procedure for it, one step further than the present paper, is described in Cersullo et al. (2019). A laser-frequency comb (LFC) was added later, that had a long period of commissioning and was offered to the community on November 15^th, 2018³. Thanks to the last upgrades of HARPS, such as the use of an LFC and octagonal fibres, the precision of the spectrograph is expected to become even better than 0.5 m s⁻¹ 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.

2 HARPS CCD block stitching

The HARPS charge-coupled device (CCD) (full detailled description in ESO 2001) is a 4K × 4K CCD mosaic made of two 2K × 4K E2V CCDs type 44–82 with a nominal pixel size of 15 μm. As shown in Fig. 1, these individual 2K × 4K CCDs are themselves made of two rows of eight 1024 × 512 pixels subparts hereafter referred to as “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 may 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 variation⁴, hence they will not be accounted for in our work.

The average dispersion of the HARPS spectrograph is 820 m s⁻¹ px⁻¹ (2.7 ×10⁻⁶ in relative) at the middle of the orders, that is 0.015 Å px⁻¹ for a wavelength of 5500 Å. Earth’s yearly revolution around the solar system barycentre (SSB) has a peak-to-peak velocity amplitude of 59.58 km s⁻¹, i.e. 2.0 × 10⁻⁴ times the speed of light in vacuum. This is equivalent to a wavelength shift of 1.1 Å at λ = 5500 Å, that is, to a position shift of 73 px on HARPS’s CCD, or about one seventh of the width of a block. The 72 HARPS spectral orders are along the 512-pixel-side direction, so each of them crosses seven 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 into the final data of a number of stars observed with this instrument. This effect has been known for along time (see e.g. Lovis & Pepe 2007, hereafter referred to as LP07, and Dumusque 2012) and is studied by several other authors (see Wilken et al. 2010; Molaro et al. 2013b; Dumusque et al. 2015; Bauer et al. 2015).

In orderto 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 an 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, we note 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. 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). Despite this fact, that solution 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, to obtain 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.

Fig. 1 Schematic representation of HARPS CCDs, highlighting the stitching pattern. The row and stitching (gap) numbers are those used in Table 1. 512 × 1024 is the size of the blocks in number of (physical) pixels.

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 m s⁻¹. The origin of this high scatter is unexplained, and knowledge of where it comes from could allow us 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 cm s⁻¹-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 Å, that had been obtained with the McMath-Pierce 1 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⁻¹⁰ relative precision aimed at by sub-m s⁻¹ RV spectrographs. LP07 took advantage of HARPS’s stability 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 ionised and doubly ionised thorium (Th I, Th II, and Th III, respectively) from 250 to 5500 nm (40 000 to 1800 cm⁻¹). They experimentally measured more than 1600 thorium lines from the near ultraviolet to the near infrared with a resolution up to 10⁷ 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 have therefore used 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 optimised 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.

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 was visible in the form of 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 is considered as two different instruments. The wavelength calibration improvement presented in this paper was therefore carried out twice independently. In the data analysis we also considered the pre-2015 and post-2015 data as two separate datasets.

5 Calibration data sets

Our work uses and is based on several different calibration data sets, each for a specific task. We describe them below for a better overview.

5.1 LED flat-field frames

The first type of frames that was used is light-emitting-diode (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 Sects. 6.1 and 6.2 respectively). The interest of using LED flat-field frames here is that they provide information that is independentfrom the wavelength calibration frames (that use a spectral source), so a consistency check can be made after correction.

As presented in Fig. 2, on these LED frames, unsurprisingly we 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). This assumes that all the pixels have the same size, and flat-field frames are, in fact, usually usedto 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 pixel size. Nevertheless, a regular pattern located at the boundaries between the blocks would be a clear indication for the latter: in fact, every 512 px we find 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 white noise. These columns are actually four to eight 1024-pixel-high sub-columns (four times 1024 px 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 Sect. 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.

Fig. 2 LED flat-field image of a part of a CCD highlighting gaps at block boundaries.

5.2 Laser-frequency-comb calibration frames

A set of 82 LFC calibration frames taken with HARPS between 2015-04-17T02:57:32.992 and 2015-04-17T04:20:07.352, not long before the May 2015 upgrade (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 Sect. 5.3) is that the absolute wavelength of each line of the laser frequency comb is known with extreme precision. Indeed, only two parameters are needed to describe the wavelengths of the whole set of laser lines: ν₀, the anchor frequency of our laser comb, and Δν, the mode spacing. Then, the formula that gives the wavelength of a given line is ν_n = ν₀ + n × Δν, for every integer n. In our case, ν₀ = 288.059800 THz and Δν = 18.0 GHz (Lo Curto, priv. comm.). 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).

5.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 the HARPS reference calibration data before and after its upgrade in May 2015 (Lo Curto et al. 2015), respectively. The first set is composed of 64 frames taken in December 2004, all with thorium-argon on both fibres, and therefore appropriate for the calibration of both fibres. The second set consists of 120 frames taken in June 2015, distributed as follows: 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 appropriate only for the calibration of fibre A. Finally, the third set contains ten frames taken between April 13th and April 23rd, 2015, used to compare the results with the ones obtained with the laser frames, that are from April 17th, 2015 (see Sect. 5.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 obtained lists of lines for the two fibres, both before and after the upgrade of May 2015. Before the upgrade, we counted a list of 10 227 lines for fibre A and a list of 10 006 lines for fibre B, most of these lines 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).

6 Gap properties determination

6.1 Positions of the gaps

To determine the positions of the gaps, we used the LED flat-field frames presented in Sect. 5.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 Fig. 2. We note that, as already stated by Molaro et al. (2013b), the positions of the large brightness anomalies are offset from the exact multiples of 512. Molaro et al. (2013b) noted an offset of four pixels between the physical positions of the gaps on the CCD (at n_j = 512j px) and the positions observed on the flat fields (at n_j = (512j + 4) px), and that one or two pixels are involved at each time. Here we find three and/or four pixels of offset, given that several of the anomalies are two pixels wide. 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).

6.2 Gap sizes

Obviously, even if the abovementioned flat-field frames already have all individually high signal-to-noise ratios (S/N), 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 Sect. 5.1 into one “master” flat-field frame.

The first hypothesis that was made is that the gaps between two given (512 × 1024)-pixel blocks have a constant size. This may not be strictly true, although 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 assumed that, by default, the pixels have a nominal size L_nom, equal to 15 μm, in both directions. The only size anomaly that we take into account is the one that we want to characterise, in other words, the gap size anomaly along the spectral orders. This assumption comes from the aforementioned 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 Fig. 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 blocks of 1024 lines of height as one single line, significantly increasing the S/N. One should notice that, assuming a random distribution of the pixel characteristics (size, sensitivity, ...), this line stacking allows us 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 seven 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 $$\frac{L_{\text{gap}}}{L_{\text{nom}}}=\frac{L_{\text{est}}}{L_{\text{nom}}}-1=\frac{F_{\text{meas}}}{F_{\text{exp}}}-1,$$(1)

where L_est is the estimated actual size of the affected pixel, F_meas is the flux actually measured on this pixel, and F_exp is the value that we would expect if the pixel had a nominal size. We note that, with this convention, the negative L_gap correspond to pixels smaller than nominal, that is, to pixels darker than normal, while the positive L_gap correspond to pixels larger than nominal, that is, 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[∪]n + 2, n + 10].

We therefore obtained 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 mpx and +15.27 mpx, corresponding to −28.3 m s⁻¹ and +12.5 m s⁻¹. Over two pixels corresponding to the same gap, the maximum values are −62.00 mpx and +25.76 mpx, i.e. −50.8 m s⁻¹ and +21.1 m s⁻¹.

We can compare these values with those obtained by Bauer et al. (2015). All except 4 are within 6 mpx (≈ 5 m s⁻¹) of the values found by Bauer et al. (2015); most of them are within 3 mpx ( ≈2.5 m s⁻¹). The average and the dispersion of the differences (Bauer et al. 2015 minus our values) are respectively −0.17 mpx (≈ − 0.14 m s⁻¹) and 6.23 mpx (≈ 5 m s⁻¹).

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, in other words, 72 lines of 4096 values. This map is used later in this work, when we compute the new wavelength solutions. The process we used to build this map is as follows. For the orders that belong to only one of the four rows of CCD blocks (which 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 used an averaged gap size weighted by the number of pixels extracted in each row at the given position on the order.

7 Wavelengths redetermination

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, for this paper, we took the gap sizes into account in the equations and that we used the new thorium wavelengths of RNS14.

7.1 Mathematical formula

Our 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 $${\displaystyle \sum _{k=0}^P{a}_k}{X}_i^k={\lambda }_i,$$(2)

where:

• {a_k}_k∈[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 i on the CCD along the order (x-axis);

• λ_i

  (known) wavelength of line i from RNS14.

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: $$X_i=x_i+{\displaystyle \sum _{j=1}^G{\delta }_{i,j}}{g}_j,$$(3)

where:

• 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;

• δ_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 i and gap j, ${\delta }_{i,j}=\{\begin{array}{ll}1\hfill & \text{if\hspace{0.17em}\hspace{0.17em}}{n}_i>n_j\hfill \\ 0\hfill & \text{else};\hfill \end{array}$

• g_j

  size of gap number j.

The ${\left\{n_j\right\}}_{j\in [1,G]}$ and the ${\left\{g_j\right\}}_{j\in [1,G]}$ are those determined respectively in Sects. 6.1 and 6.2.

When we insert Formula (3) inside Formula (2), we get $${\displaystyle \sum _{k=0}^P{a}_k}{\left(x_i+{\displaystyle \sum _{j=1}^G{\delta }_{i,j}}{g}_j\right)}^k={\lambda }_i.$$(4)

The only unknowns are the ${\left\{a_k\right\}}_{k\in [0,P]}$, and Eq. (4) is then linear with respect to the unknowns. For simplicity, from here on we use the term X_i (= $x_i+{\displaystyle \sum _{j=1}^G{\delta }_{i,j}}{g}_j)$.

Practically, Eq. (4) will be different depending on the type of line that we are dealing with. We have four 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).

• 1.

  Duplicate thorium lines: $${\displaystyle \sum _k}{a}_k{X}_i^k-\Delta {\lambda }_i={\lambda }_{i,\text{\hspace{0.17em}RNS14}}.$$(5)

• 2.

  Unique thorium lines: $${\displaystyle \sum _k}{a}_k{X}_i^k={\lambda }_{i,\text{\hspace{0.17em}RNS14}}.$$(6)

• 3.

  Duplicate non-thorium lines: $${\displaystyle \sum _k}{a}_k{X}_i^k-{\lambda }_i=0.$$(7)

• 4.

  Unique non-thorium lines: no equation

In these equations, the Δλ_i is the unknown systematic offset of line i (for thorium lines), λ_{i, RNS14} is the RNS14 wavelength of line i (for thorium lines) and λ_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 Δλ_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 way also strengthen the reliability of the wavelength solutions towards the order edges, where the wavelength solutions are usually more poorly constrained.

Fig. 3 Residuals of the duplicate thorium lines. The black line is the line of equation y = x.

7.2 Thorium-line systematics

The individual systematics Δλ_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 theleft-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 Δλ_i) are on average equal to zero; second, the Δλ_i of the equations corresponding to a same duplicate line should all be considered to be equal. In addition, it means that the Δλ_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 Δλ_i. These global properties of the Δλ_i translate into these additional equations (one for each Δλ_i, i.e. one for each duplicate thorium line): $$\Delta {\lambda }_i=0.$$(8)

This has to be understood in a least-squares sense: these equations will ensure that, in the global solution, the properly normalised (see below) Δλ_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. Sect. 7.3).

7.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 equations of type (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 px). 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 see in Sect. 8.1.2 and Fig. 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 Δλ_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. We note that, here, the flux are the flux measured by HARPS, corrected from the blaze and averaged 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 ofboth quantities. This is shown in Fig. 4, where we used ten 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 Δλ_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: $$\log (\sigma [\AA ])=-0.2488\times \log (I [{\text{e}}^{-}])-\mathrm{1.8311.}$$(9)

i.e. $$\sigma [\AA ]=0.01475\times {(I [{\text{e}}^{-}])}^{-0.2488}.$$(10)

For a line at 5500 Å, this relation corresponds to an uncertainty on the systematic offset equivalent to 25 m s⁻¹ for a flux of one million electrons and to about 150 m s⁻¹ 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 m s⁻¹. 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 the picture 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.

Fig. 4 Dispersion of the thorium-line systematic offsets Δλi (in Å) vs. logarithm of the line intensity (in electrons).

7.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: $$\lambda =\frac{K}{m},$$(11)

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 use a third-degree polynomial: $$\lambda ={\sum }_c{b}_c{\left(\frac{1}{m}\right)}^c.$$(12)

For a given x and a given order, we know that this point must satisfy both the wavelength solution of its order (Eq. (2)) and the Littrow equation of its position (Eq. (12)): $${\sum }_k{a}_k{X}^k-{\sum }_c{b}_c{\left(\frac{1}{m}\right)}^c=0.$$(13)

In orderto 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 }_{\text{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: $$N\times \frac{1}{{\sigma }_{\text{L}}^2}={\sum }_i\frac{1}{{\sigma }_i^2}.$$(14)

The corresponding uncertainty attributed to these equations is therefore, $${\sigma }_{\text{L}}=\sqrt{\frac{N}{{\displaystyle \sum _i(1/{\sigma }_i^2)}}}.$$(15)

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

7.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 Sect. 5.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 σ away from the solution of the first run. We used the solution of this second run 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 value 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 yr. This new list is specific to HARPS and contains the updated wavelengths that include the line systematics Δλ_i. The typical dispersion of the residuals around the new wavelength solutions is 10 m s⁻¹ for the pre-upgrade data and 15 m s⁻¹ for the post-upgrade data. While the pre-upgrade solutions show similar dispersions compared to the old reduction, we observe a significant improvement in the post-upgrade solutions (from 30 to 15 m s⁻¹).

8 Results

8.1 Test of the gap correction with the laser frequency comb

8.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 λ = P(x), where λ 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. 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 m s⁻¹, larger than the RV-converted size of the gaps. The set of 82 laser frames mentioned above (Sect. 5.2) was therefore used to perform this test.

As mentioned above, the laser comb does not currently cover the whole CCD, but only orders 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 row 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.

8.1.2 Results

The position of the lines were first fitted in pixels, meaning 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 Fig. 5, by doing this correction, we do not see any noticeable jump, which indicates that the estimated sizes of the gaps are correct within the residual dispersion of about 7 m s⁻¹.

This low residual dispersion of the laser solutions, about 7 m s⁻¹, i.e. about 9 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 m s⁻¹) observed with thorium frames. It reinforces the suggestions that the large scatter seen on thorium solutions is mainly due to intrinsic features of the thorium spectrum.

8.2 Comparison of the thorium and laser solutions

We further compared the laser and thorium wavelength solutions. In Fig. 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.

We do, however, notice that the new wavelength solution has an average offset of about 43 m s⁻¹ with respect to the laser solution, while the standard thorium solution has an offset of about 49 m s⁻¹. The offset seems to be larger for longer wavelengths. As we explained in Sect. 5.2, the theoretical wavelengths of the LFC respect the formula ν_n = ν₀ + n × Δν, where Δν (18 GHz in our case) is known at the Hz level and ν₀ was set at 288.0598 THz (LFC team, priv. comm.). We note that a fixed error Δν₀ on ν₀ would translate into an error in RV proportionalto the wavelength: Δv = Δν₀ × λ. An error Δν₀ = 100 MHz would translate into an offset Δv = 40, 50, and 60 ms⁻¹ for a wavelength λ = 4000, 5000, and 6000 Å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, ν₀ was varied by a few hundred MHz over the course of the commissioning. A mistake in the reported ν₀ for the specific frames we used 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 ν₀ 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 ν₀ 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 ν₀ 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 note that, if we compare the two thorium wavelength solutions, we have a mean offset of $c\times \frac{{\lambda }_{\text{new}}-{\lambda }_{\text{old}}}{{\lambda }_{\text{laser}}}\approx -6$ m s⁻¹. This is fully compatible with the correction factor κ = 0.22 × 10⁻⁷ applied by RNS14 on PE83’s wavenumbers, that corresponds to an offset of $c\times \frac{{\lambda }_{\text{RNS14}}-{\lambda }_{\text{PE83}}}{{\lambda }_{\text{RNS14}}}=-6.6$ m s⁻¹. This has to be expected since the re-computation of the wavelengths by LP07 and by us are fit, for the absolute values, on PE83 and RNS14 respectively.

Fig. 5 Comparison of the residuals of laser frequency comb wavelength solutions, with or without correction of the gaps, for orders 25, 42 and 57.

Fig. 6 Difference, for the same orders as Fig. 5, of the wavelength solution from the laser frequency comb and that from the thorium, i.e c0 × (λTh − λLFC)∕λLFC, where c0 is the speed of light in vacuum, for both the standard (green) and new (red) thorium solutions. The new thorium wavelength solutions match the laser frequency comb solutions much more closely. Especially, we can see that there are “jumps” on the green curves, while the red curves are continuous. We note that the laser solutions are fifth-order polynomials while the thorium solutions are third-order polynomials.

Fig. 7 Cumulative histogram of the rms of the HARPS RVs of all the objects observed at least10 times during the period by HARPS that have a radial velocity rms below 2 m s−1 with both the standard (blue) and the new (green) versions of the DRS.

8.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 in 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 ten times during the period by HARPS that have a radial velocity rms below 2 m s⁻¹ 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 ten 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 twenty most stable stars all have their rms improved (except one), in many cases by a significant margin. An example of this improvement is shown in Fig. 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 to 1.23 m s⁻¹. This improvement is especially due to the recomputation of the thorium wavelengths for the post-upgrade HARPS as explained in Sect. 7.

Fig. 8 HARPS RVs of HD 220507 after the change of fibres (82 points, Δt = 803 days). The rms of the data is 1.73 m s−1 with the standard version of the DRS (purple) while it reduces to 1.23 m s−1 with the new version of the DRS (orange). No significant signal is observed for this star.

8.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 yr) 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 those obtained with the standard pipeline, that did not take into account the gaps and used the wavelengths as defined in LP07. In particular, 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 one-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 Exoplanets⁵ (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 HD 10700, HD 20794, HD 69830, HD 1461, HD 154088, HD 31527, and HD 40307. For all the systems, a noise level of 0.75 m s⁻¹ 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 allowed 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 we list 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.

8.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 ~1 m s⁻¹ (Pepe et al. 2011; Tuomi et al. 2013a). 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 one-year signal becomes slightly less significant with the new version (FAP of 0.49%) compared tothe 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. (2013a) 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, which 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. (2013a) 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 demonstratethat their amplitude and significance depend on the pipeline version used and the data analysis method.

8.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 ± 0.005 d and 89.65 ± 0.10 d with the standard DRS and at 18.320 ± 0.005 d and 89.76 ± 0.12 d with the new DRS. For the amplitudes of these signals, we find respectively 0.67 ± 0.07 and 0.71 ± 0.07 m s⁻¹ for HD 20794 b, and 0.74 ± 0.09 and 0.57 ± 0.07 m s⁻¹ for HD 20794 d. While the difference is within 1σ for planet b, it is more noticeable for planet d, 22% and nearly 2σ 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 ± 0.10 (b) and 0.45 ± 0.09 (d) with the standard version of the DRS to 0.11 ± 0.10 (b) and 0.29 ± 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 thesesignals 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 inorder 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.

8.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 ≈ 3680 d, K ≈ 1.1 m s⁻¹, e ≈ 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 ± 0.0004 days, an amplitude of 3.47 ± 0.11 m s⁻¹ and an eccentricity of 0.07 ± 0.03. For the second planet, the period is 31.601 ± 0.006 days, the amplitude 2.66 ± 0.13 m s⁻¹ and the eccentricity 0.06 ± 0.05. Finally, the outermost planet is found to have a period of 200.32 ± 0.36 days, an amplitude of 1.85 ± 0.12 m s⁻¹ and an eccentricity of 0.12 ± 0.06 with the standard version of the DRS while we find a period of 200.57 ± 0.37 days, an amplitude of 1.77 ± 0.12 m s⁻¹ and an eccentricity of 0.10 ± 0.06 with the new version.

Fig. 9 Periodograms of HD 10700 without fitting any Keplerian signal or drift.

Fig. 10 Periodograms of HD 20794 after the fit of the two known planets at 18.32 and 89.7 days.

Fig. 11 Periodograms of HD 69830 after the fit of the three known planets at periods of 8.67, 31.60, and 200 days.

Fig. 12 Periodograms of HD 1461 after the fit of the two known planets at 5.77 and 13.50 days and the magnetic cycle at 3580 days.

8.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 days) around this star as well as a magnetic cycle (period of 3500 days). 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 one year, even if we notice that the amplitude is a little 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σ 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 ± 80 days, an amplitude of 3.28 ± 0.28 m s⁻¹ and an eccentricity of 0.51 ± 0.07. Concerning the innermost planet, we find a period of 5.7715 ± 0.0004 days, an amplitude of 2.26 ± 0.15 m s⁻¹ and an eccentricity of 0.11 ± 0.07. For the second planet, we get a period of 13.506 ± 0.003 days, an amplitude of 1.60 ± 0.16 m s⁻¹ and an eccentricity of 0.15 ± 0.10.

Fig. 13 Periodograms of HD 154088 after the fit of 1 planet at 18.55 days.

Fig. 14 Periodograms of HD 154088 after the fit of 1 planet at 18.55 days and a linear drift.

8.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 theHARPS 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 signal and with a Keplerian signal plus a linear drift. In both cases, all the values from the old and the new DRS are compatible within 1σ so, as before, we give only the values from the new DRS.

From the model with only a Keplerian signal, we find a period of 18.556 ± 0.009 days, nearly 2σ lower than the 18.596 ± 0.021 days in Mayor et al. (2011). The amplitude is found to be 1.68 ± 0.19 m s⁻¹, slightly lower than but well within 1σ of the 1.78 ± 0.31 m s⁻¹ of Mayor et al. (2011). We find an eccentricity of 0.10 ± 0.11, noticeably lower than the 0.38 ± 0.15 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}_{-6487}^{+\infty }$ days in this paper. Here, we 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 ± 0.07 m s⁻¹ yr⁻¹. In this case, we find a period of 18.568 ± 0.008 days for the planet, slightly larger than before but still lower than found by Mayor et al. (2011). The amplitude is found to be 1.64 ± 0.23 m s⁻¹, slightly lower than before but well within 1σ. We find an eccentricity of 0.28 ± 0.12, 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, we see that the peak around 360 days is smaller with the new DRS than with the standard one.

8.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σ, except for the amplitudes of planet d, that are nonetheless within less than 2σ. 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 ± 0.003 days, an amplitude of 2.84 ± 0.13 m s⁻¹ and an eccentricity of 0.16 ± 0.05. HD 31527 c has a period of 51.20 ± 0.03 days, an amplitude of 2.70 ± 0.14 m s⁻¹ and an eccentricity of 0.08 ± 0.05. Finally, HD 31527 d has a period of 274.2 ± 1.0 days, an amplitude of 1.77 ± 0.36 m s⁻¹ and an eccentricity of 0.67 ± 0.08. 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. We also note that, with 89 more points and nearly twice as long a time-span as in Mayor et al. (2011), we have significantly lower uncertainties on the parameters of the outermost planet.

Fig. 15 Periodograms of HD 31527 after the fit of the three known planets: 16.56, 51.2, and 272 days.

Fig. 16 Periodograms of HD 40307 after the fit of four signals at periods of 4.31, 9.62, 20.42, and 51.6 days as well as a third-degreepolynomial for the magnetic cycle.

8.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. (2013b) 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. (2013b). While they find a signal at 205 days that may correspond to Tuomi et al. (2013b)’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 isstrongly 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 remains. 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σ, so here we also give only the values from the new version. The innermost planet has a period of 4.3114 ± 0.0002 days, an amplitude of 1.84 ± 0.14 m s⁻¹ and an eccentricity of 0.12 ± 0.07. The second planet has a period of 9.6210 ± 0.0008 days, an amplitude of 2.29 ± 0.13 m s⁻¹ and an eccentricity of 0.05 ± 0.06. The third planet has a period of 20.412 ± 0.004 days, an amplitude of 2.31 ± 0.14 m s⁻¹ and an eccentricity of 0.07 ± 0.05. The fourth signal has a period of 51.68 ± 0.07 days, an amplitude of 0.78 ± 0.13 m s⁻¹ and an eccentricity of 0.19 ± 0.16. Finally, the drift parameters are 0.60 ± 0.07 m s⁻¹ yr⁻¹, − 0.13 ± 0.02 m s⁻¹ yr⁻² and − 0.02 ± 0.00 m s⁻¹ yr⁻³.

8.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 Sect. 8.3 and in Fig. 8. Concerning the fitted signals (confirmed planets and magnetic cycles), all the amplitudes are compatible within 2σ, most of them within 1σ, between the two versions of the DRS, as well as between these results and those in the literature.

We also compared the periodogram power around one year between the two versions of the DRS. Among the seven systems studied, the periodogram power at one year decreases with the new reduction in six cases. In the only case where it slightly increases (HD 31527), the power at one year remains below the significance threshold. This demonstrates that spurious instrument-related one-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 seven systems, only one (HD 1461) seems to have a significant signal at one year with the new reduction. It is possible that this particular signal is in fact of astrophysical origin.

8.5 Test on individual stellar spectral lines

To see the correction that our new wavelength solution brings at the location of the stitching 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. In the bottom panels of Fig. 17 we also show 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. 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.

Fig. 17 RVs of three spectral lines in the spectrum of HD 1461 that passes over CCD stitching boundaries. We plot those RVs as a function of phase during the year, because due to Earth orbiting around the Sun, the stellar spectrum will move on theCCD with a period of a year. This implies that the spectral lines selected here will pass over CCD stitching boundaries, shown in these plots with the red vertical dashed lines, always at the same phase. The top panels show the RV derived using the standard DRS and the middle panels the RV derived with the correction described in this paper. We show in the bottom panels the difference between those two RV sets. The blue horizontal lines represent the median RV before and after the spectral line passes over the CCD stitching boundary. As we can see, the significant offset between the DRS RVs before and after the CCD stitching boundary disappears in the RVs derived with the correction described in this paper, proving that we consider correctly the different pixel sizes at the stitching positions.

9 Conclusions

The aim of this work was to improve the wavelength calibration of the HARPS spectrograph to push further its planet-detection capabilities.

• 1.

  We estimated 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 px. The validity of this correction was verified with laser frequency comb exposures.

• 2.

  We reprocessed 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 block stitching boundaries, while the old solutions show significant jumps.

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

  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.

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

• 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 detail in Cersullo et al. (2019).

Acknowledgements

We thank the anonymous referee for his/her remarks that greatly helped to 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.

References

All Tables

All Figures

	Fig. 1 Schematic representation of HARPS CCDs, highlighting the stitching pattern. The row and stitching (gap) numbers are those used in Table 1. 512 × 1024 is the size of the blocks in number of (physical) pixels.
In the text

	Fig. 2 LED flat-field image of a part of a CCD highlighting gaps at block boundaries.
In the text

	Fig. 3 Residuals of the duplicate thorium lines. The black line is the line of equation y = x.
In the text

	Fig. 4 Dispersion of the thorium-line systematic offsets Δλi (in Å) vs. logarithm of the line intensity (in electrons).
In the text

	Fig. 5 Comparison of the residuals of laser frequency comb wavelength solutions, with or without correction of the gaps, for orders 25, 42 and 57.
In the text

Fig. 6 Difference, for the same orders as Fig. 5, of the wavelength solution from the laser frequency comb and that from the thorium, i.e c0 × (λTh − λLFC)∕λLFC, where c0 is the speed of light in vacuum, for both the standard (green) and new (red) thorium solutions. The new thorium wavelength solutions match the laser frequency comb solutions much more closely. Especially, we can see that there are “jumps” on the green curves, while the red curves are continuous. We note that the laser solutions are fifth-order polynomials while the thorium solutions are third-order polynomials.

In the text

	Fig. 7 Cumulative histogram of the rms of the HARPS RVs of all the objects observed at least10 times during the period by HARPS that have a radial velocity rms below 2 m s−1 with both the standard (blue) and the new (green) versions of the DRS.
In the text

	Fig. 8 HARPS RVs of HD 220507 after the change of fibres (82 points, Δt = 803 days). The rms of the data is 1.73 m s−1 with the standard version of the DRS (purple) while it reduces to 1.23 m s−1 with the new version of the DRS (orange). No significant signal is observed for this star.
In the text

	Fig. 9 Periodograms of HD 10700 without fitting any Keplerian signal or drift.
In the text

	Fig. 10 Periodograms of HD 20794 after the fit of the two known planets at 18.32 and 89.7 days.
In the text

	Fig. 11 Periodograms of HD 69830 after the fit of the three known planets at periods of 8.67, 31.60, and 200 days.
In the text

	Fig. 12 Periodograms of HD 1461 after the fit of the two known planets at 5.77 and 13.50 days and the magnetic cycle at 3580 days.
In the text

	Fig. 13 Periodograms of HD 154088 after the fit of 1 planet at 18.55 days.
In the text

	Fig. 14 Periodograms of HD 154088 after the fit of 1 planet at 18.55 days and a linear drift.
In the text

	Fig. 15 Periodograms of HD 31527 after the fit of the three known planets: 16.56, 51.2, and 272 days.
In the text

	Fig. 16 Periodograms of HD 40307 after the fit of four signals at periods of 4.31, 9.62, 20.42, and 51.6 days as well as a third-degreepolynomial for the magnetic cycle.
In the text

Fig. 17 RVs of three spectral lines in the spectrum of HD 1461 that passes over CCD stitching boundaries. We plot those RVs as a function of phase during the year, because due to Earth orbiting around the Sun, the stellar spectrum will move on theCCD with a period of a year. This implies that the spectral lines selected here will pass over CCD stitching boundaries, shown in these plots with the red vertical dashed lines, always at the same phase. The top panels show the RV derived using the standard DRS and the middle panels the RV derived with the correction described in this paper. We show in the bottom panels the difference between those two RV sets. The blue horizontal lines represent the median RV before and after the spectral line passes over the CCD stitching boundary. As we can see, the significant offset between the DRS RVs before and after the CCD stitching boundary disappears in the RVs derived with the correction described in this paper, proving that we consider correctly the different pixel sizes at the stitching positions.

In the text