© ESO 2021

1. Introduction

The radial-velocity (RV) method is an efficient technique to estimate the mass of exoplanets. If the planets are also transiting, their radius and exact mass can be obtained, thus revealing the exoplanet’s density, leading to a planetary interior and atmospheric characterisation. The RV technique is also currently the only method for which available instruments reach the sensitivity needed to detect low-mass, non-transiting exoplanets orbiting stars within 40 parsecs from the Sun. Those planets will be prime targets for future direct imaging missions (e.g., National Academies of Sciences, Engineering, Medicine 2018), such as LUVOIR¹ and HabEx².

However, when looking for Earth-mass planets orbiting in the habitable zone of solar-type stars, the RV technique is limited by spurious radial-velocity signals induced by the host star at the m s⁻¹ level (e.g., Fischer et al. 2016; Plavchan et al. 2018; Blackwood et al. 2020). It is therefore crucial to obtain a better understanding of stellar signals in order to mitigate them and thus reveal tiny planetary signatures with amplitudes as small as 10 cm s⁻¹.

Due to stellar signal timescales ranging from minutes for stellar oscillations to years for magnetic cycles, standard radial-velocity observations obtained either from transit follow-up observations or blind search surveys very often lack the sampling and time baselines necessary to characterise the different types of stellar signals. The lack of good data makes it extremely difficult to understand stellar signals in enough detail and, therefore, to find efficient techniques to mitigate them. The ideal solution to make progress is to continuously observe a target for which astrophysical parameters and activity behaviour as a function of time are perfectly known, with the same instrument employed for exoplanet characterisation. This reflection led to the development of the HARPS-N low-cost solar telescope (LCST, Dumusque et al. 2015). Since then, several solar feeds have been installed at different sub- m s⁻¹ RV precision facilities and others are planned for the future (HARPS, EXPRES, NEID, MAROON-X, priv. comm. with the different instrument teams). We note that a solar feed is also present on the PEPSI spectrograph (Strassmeier et al. 2015) and on the NIR spectrograph GIANOB (Claudi et al. 2018); however, the obtained data are not dedicated to extreme RV precision.

Since July 2015, HARPS-N has been observing the Sun for several hours nearly every day, with a cadence of 5 min, using its solar feed (Phillips et al. 2016). Since then, several papers have analysed the obtained data and investigated solar radial-velocity signals and related activity-index indicators (Collier Cameron et al. 2019; Milbourne et al. 2019; Maldonado et al. 2019; Miklos et al. 2020; Langellier et al. 2020)

The excellent quality of the solar RVs published in Collier Cameron et al. (2019), with instrumental systematics below the m s⁻¹ level, and the exceptional cadence of these data over several years, revealed instrumental systematics that needed to be understood and corrected for. The first step in this direction was to adapt the publicly available ESPRESSO Data Reduction Software (DRS, Pepe et al. 2021) to HARPS-N so that the solar data could be reduced with an up-to-date DRS. In this paper, we refer to the current HARPS-N DRS (version 3.7) as the old DRS and the ESPRESSO pipeline (version 2.2.3) adapted to HARPS-N as the new HARPS-N DRS. Beyond using the ESPRESSO DRS to extract RV information from raw spectra, we also developed new recipes to specifically account for HARPS-N systematics. We implemented a new procedure to derive wavelength solutions, optimised the selection of thorium lines to improve RV stability in the long term, used master flat fields to prevent photon noise being limited by them when analysing solar data, used the calibrations that are the closest in time, and corrected for ghost’s contamination in the calcium activity index time series.

The article is structured as follow. In Sect. 2 we discuss the optimisation of the ESPRESSO DRS to HARPS-N and the new recipes implemented to solve for known systematics. In Sect. 3, we discuss the HARPS-N solar data set used in this paper, we compare the solar RVs reduced with the old and new DRS, and we describe the data that are released alongside this paper. Finally, we discuss the results found and conclude in Sect. 4.

2. Deriving precise HARPS-N radial-velocity measurements

2.1. The ESPRESSO DRS

To improve the precision of the HARPS-N RVs, we first adapted the publicly available ESPRESSO DRS to work with HARPS-N. Although many of the data reduction steps and algorithms implemented in the ESPRESSO DRS have origins in the HARPS-N DRS, several steps were significantly re-designed and improved with respect to HARPS-N to meet the more challenging requirements of ESPRESSO. A comparison of the old HARPS-N and new ESPRESSO DRS is beyond the scope of this paper, and interested readers are encouraged to read Pepe et al. (2021) for a general description and, in the near future, the paper specifically describing the ESPRESSO DRS.

2.2. Improved wavelength solution

In the old HARPS-N DRS, a wavelength solution for the science fibre is performed every day by illuminating the corresponding fibre with a thorium-argon lamp. A list of thorium lines, 3692 in total, with their physical position on the detector and their laboratory wavelength is used to detect and measure the centroid of each thorium line in the extracted calibration spectrum. Readers interested in how those lines have been selected should refer to Lovis & Pepe (2007). Using the line centroids and corresponding laboratory wavelengths, a wavelength for each pixel is derived by interpolation. For each order individually, this interpolation is performed using a third order polynomial. The old HARPS-N DRS, similar to the new one, extracts 69 orders from the raw images; thus, wavelength solutions are derived using a total of 4 × 69 = 276 free parameters.

We note that to ease the comparison with the DRS data products, we use the order numbering of the DRS throughout this paper. Order 0 thus corresponds to the first row in the echelle-order spectra derived from the DRS, order 68 to the last one. In terms of physical orders, order 0 of the DRS corresponds to physical order 158 (average wavelength 6873 Å) and order 68 corresponds to physical order 90 (average wavelength 6872 Å).

The flux throughout an echelle order is not constant due to the blaze response, and it is very low on both sides of the order. Therefore, thorium lines in the centre of an order are much more constraining for the signal-to-noise (S/N) weighted polynomial fit used to derive the wavelength solution than thorium lines on both sides of the order. This gives rise to instability in the fitted polynomials and therefore in the wavelength solution from calibration to calibration, which are generally performed every day. This problem has already been discussed in several papers (Dumusque 2018; Cersullo et al. 2019; Bauer et al. 2015) and in the top row of Fig. 1 we show an example of how the wavelength solution for a few orders can differ by tens of meters-per-second from one calibration to the next. Here, we show several wavelength solutions spread over 200 days, and we observe that these strong variations can be seen from one calibration to the next, even when taken only a few minutes apart. This instability in wavelength solutions translates into day-to-day offsets in RV measurements with an estimated root-mean-square (rms) of about 0.4 m s⁻¹ (Dumusque 2018).

Fig. 1. RV difference between wavelength solutions taken on different days when derived using different algorithms. The RV difference was obtained by comparing all the wavelength solutions with the wavelength solution of reference used in the drift map method, obtained on 2017 November 11. Each column corresponds to a different echelle order, highlighting the behaviour of the wavelength solutions from the blue (order 3) to the red (order 60). Top: RV difference when using only the thorium lines (TH). This is what is done in the old HARPS-N DRS. Middle: RV difference when using thorium lines and the Fabry-Pérot spectrum (TH–FP). Bottom: RV difference when using the drift map method presented in this paper, which means using a reference wavelength solution, and shifting it in velocity after measuring the drift of the thorium lines.

To solve for this instability issue in the wavelength solution, the ESPRESSO DRS uses a combination of a thorium calibration on the science fibre, followed by a Fabry-Pérot calibration on the same fibre to stabilise the wavelength solution. The thorium lines give the absolute wavelength scale, while the Fabry-Pérot’s very rich spectrum, with each emission peak being separated by the same distance due to the interferometric nature of the emission, allows us to better constrain the polynomial fit. Readers interested in how to implement such a technique are directed to Cersullo et al. (2019) and Bauer et al. (2015). However, for HARPS-N, such a set of calibrations is not present during the entire lifetime of the instrument. As we wanted to adopt a single wavelength solution derivation method for the spectrograph, we decided to implement another strategy that we call the drift map method. The idea is to use a reference wavelength solution, obtained at a time when consecutive thorium and Fabry-Pérot calibrations on the science fibre are available, and then measure the drift of the thorium lines between each thorium calibration and the reference thorium calibration. A two-dimensional polynomial fit is then performed on the measured drift of the thorium lines and a wavelength solution is obtained for each thorium calibration by shifting the reference wavelength solution accordingly. A detailed description on how to derive a wavelength solution using the drift map method can be found in Appendix B.

To test how the drift map method compares to the thorium-only (old HARPS-N DRS) and combined thorium–Fabry-Pérot techniques when deriving wavelength solutions, we selected 29 thorium calibrations on the science fibre, followed by Fabry-Pérot calibrations on the same fibre (from 2017 March 25 to 2017 November 06). We then derived the wavelength solutions for those calibrations using the thorium–Fabry-Pérot technique and the drift map method, and we compared the obtained products to the wavelength solutions derived using the old HARPS-N DRS. In Fig. 1, we show the RV drift, in meters-per-second between each of these wavelength solutions and the reference wavelength solution used in the drift-map method, obtained on 2017 November 11. From top to bottom, we show the RV drift for the wavelength solution derived using the old HARPS-N DRS, the thorium–Fabry-Pérot combined solution, and the drift map method, respectively. Each column in the figure corresponds to different orders from the blue to the red, including echelle orders 3, 7, 51, and 60. We note that during those 8 months, the instrument drifted by about +15 m s⁻¹, which induces an increasing RV offset when going from the 2017 March 25 to the 2017 November 06 wavelength solutions (purple to yellow coloured lines in Fig. 1). However, we measured this instrumental drift on the drift map wavelength solutions and corrected it for all wavelength solutions. To obtain a mean drift value for each drift map wavelength solution, we performed a weighted average taking into account the blaze response of the instrument and the change in S/N across the echelle orders.

Looking at Fig. 1, it is clear that the thorium-only wavelength solutions show more variations from one wavelength solution to the next than the thorium–Fabry-Pérot and drift map wavelength solutions for orders 51 and 60. For blue orders 3 and 7, the thorium–Fabry-Pérot solutions show similar variations as the thorium-only solutions, and only the drift map method gives more stable results. This is due to the fact that the flux on the HARPS-N Fabry-Pérot is rather low in the ten first orders of the spectrograph. Fabry-Pérot emission peaks on the borders of these blue orders are difficult to detect due to low S/N, and therefore we needed to add the information provided by the thorium lines to improve the wavelength solutions in the blue. Without adding the information from the thorium lines, the thorium–Fabry-Pérot wavelength solutions in the blue part of HARPS-N are even worse than the thorium-only solution. We note that the HARPS-N Fabry-Pérot was upgraded after 2017 November 11 and the Fabry-Pérot spectrum now has much more flux in the blue part of the spectrograph, as is also the case for ESPRESSO. For all the orders, from the extreme blue to the extreme red, we see that the wavelength solutions derived with the drift map method are the most stable ones. We note that although we see that the drift map wavelength solutions present the smallest variations, this is biased by the fact that rather than fitting 276 free parameters to derive a wavelength solution, our drift map model only includes 15 parameters (see Appendix B). In the absence of an absolute wavelength solution derived using a laser frequency comb, it is impossible to say which of the thorium–Fabry-Pérot or drift map wavelength solutions are the more precise for echelle orders greater than ∼10. Since the drift map method gives wavelength solutions that are much more stable than the old DRS for all the echelle orders, we decided to implement this technique in the new HARPS-N DRS.

Wavelength solutions set the RV zero-point of the spectrograph each time a new calibration is obtained, which is every day on average. Thus, their precision directly impacts the stability of the measured RVs from one day to the next. To measure the gain in RV stability brought by the new drift map wavelength solutions compared to what was done in the old HARPS-N DRS, we measured the mean RV drift of all wavelength solutions relative to the reference wavelength solution used (set for 2017 November 11, see above). As in the analysis performed above, the mean RV drift is measured using a weighted average taking into account the blaze response of the instrument and the varying S/N across echelle orders. In the top panel of Fig. 2, we show the result for the old and new wavelength solutions. We clearly see that the old and new DRS give a similar HARPS-N drift over three years of a hundred meters-per-second. We also observe sudden jumps at each instrument interventions, either corresponding to a scheduled detector warm-up or to a power failure. In both cases, the observed offset likely occurs because after warming up and cooling down the detector, it does not come back exactly at the same previous position (see Sect. 2.6 and Fig. 5 for the rationale behind periodically warming up the detector).

Fig. 2. RV stability of the wavelength solutions derived with the old (blue) and new (orange) DRS. Top: RV drift of HARPS-N over time, measured by computing the RV drift between all wavelength solutions and the wavelength solution of reference (set for 2017 November 11, see text). We see a general drift over time, interrupted by significant jumps at each warm-up of the detector and power failures (dashed and dotted vertical lines). Bottom: distribution of the RV offsets measured between consecutive wavelength solutions separated by less than 1.5 days, which removes the large offsets observed at each instrument intervention. The MAD of these day-to-day RV offsets is 0.74 and 0.49 m s−1 for the old and new DRS, respectively.

Over the long term, several effects can influence the RV stability of the spectrograph, including the ageing of the thorium-argon lamp and the change of the focus of the instrument due to ageing of the optics. In the absence of an absolute calibrator as a laser frequency comb, we can only measure the precision in the wavelength solution over short timescales, considering that the instrument did not drift. To do so, we measured the RV offset between consecutive wavelength solutions that are separated by less than 1.5 days. This selection in time separation allows one to reject the large RV offsets measured between wavelength solutions taken before and after an instrument intervention. In the bottom panel of Fig. 2, we show the distribution of these measured day-to-day RV offsets. As we can see, the new DRS gives a distribution that is more peaked around zero and that presents less outliers. The median absolute deviation (MAD³) of the two distributions is 0.74 and 0.49 m s⁻¹ for the old and new wavelength solutions, respectively. Consequently, by using the drift map method to derive wavelength solutions, we improved the stability of the RV measurements by a non-negligible $\sqrt{0.{74}^2-0.{49}^2}=0.55$ m s⁻¹. We note that in Dumusque (2018, see Appendix C), a similar comparison using HARPS stellar data showed that the drift map wavelength solutions improved the RV precision by 0.4 m s⁻¹, a value slightly lower than what is found here for HARPS-N.

2.3. New set of thorium lines for wavelength solutions

After fitting for the position of the thorium lines, the old and new DRS check if lines are saturated and compare the fitted positions with reference positions. Lines that show a position significantly different from the reference are rejected to prevent biasing the final wavelength solution. With ageing of the thorium lamp over time, more and more lines are rejected due to lines that start to saturate and lines for which blends change in intensity or in position⁴, which induce a line-profile asymmetry and therefore a spurious RV drift. In the left panel of Fig. 3, we can see that out of a total of 3332 thorium lines used by the old HARPS-N DRS on 2014 August 16 (BJD = 2456886), only 2794 (16% less) were still used 5.5 years later⁵ on 2019 January 16 (BJD = 2458500). In addition, we can see in Appendix C that the old HARPS-N DRS was using thorium lines for which the measured RV is strongly correlated with the flux variation of the thorium-argon lamp, thus pointing towards spurious RV shifts. Finally, the old HARPS-N DRS fitted a combination of several Gaussian profiles on main thorium lines to account for blends, as described in Lovis & Pepe (2007), while the new HARPS-N DRS only fits a single Gaussian profile for each thorium spectral line. The fact that the number of thorium lines used changes over time and that the new DRS does not fit for blends anymore can have a significant impact on the long-term RV stability. We therefore decided to perform a tailored selection of thorium lines for the new HARPS-N DRS. The description of this selection is rather technical and can be found in Appendix C.

Fig. 3. Left: number of thorium lines used over time to perform the wavelength solution in the old and new HARPS-N DRS. Here, we show the results for all thorium–thorium and thorium–Fabry-Pérot calibrations. We see that there is a significant variation in the number of lines used in the old HARPS-N DRS, 16%, while in the new HARPS-N DRS, we nearly always used 2041 thorium lines. Right: RV difference between the wavelength solutions obtained with the old HARPS-N DRS and the new one, which uses the drift map method with the new set of thorium lines.

After carefully selecting thorium lines, we were left with 2041 lines that could be used to derive wavelength solutions based on the drift map method. Those thorium lines are nearly always used throughout the lifetime of HARPS-N (99% of the time, see Appendix C). We note that the entire study performed on the drift map wavelength solutions in Sect. 2.2 uses this new set of thorium lines.

In the right panel of Fig. 3, we plotted the RV difference between the wavelength solutions obtained with the old HARPS-N DRS and the new one, which uses the drift map method with the new set of thorium lines. Overall, we see an offset of 4 m s⁻¹ between the wavelengths solutions, which impacts the accuracy⁶ of the derived RVs, but not their precision, which is generally what we want. At the second order level, we see that the RV difference between the wavelength solutions is slightly correlated (Pearson correlation coefficient of 0.39) with the thorium flux ratio (see top panel of Fig. C.1), which measure the total flux emitted by the thorium-argon lamp at a given time relative to a reference time. As the thorium flux ratio significantly increases after BJD = 2458070, the difference in RV is also affected and presents a slope of ∼0.7 m s⁻¹ yr⁻¹. This tells us that the old or the new DRS seems more sensitive to flux variations of the thorium-argon lamp, however, without telling us which one it is. This can only be done by looking at the long-term monitoring of standard stars, which we do for the Sun in Sect. 3.2.

2.4. Master flat fielding

In standard HARPS-N calibrations, five consecutive frames with a tungsten lamp illuminating both science and reference fibres are used to derive the flat field for each echelle order and the corresponding blaze. For most of the night-time observations, even when a few spectra are combined together, those flat fields have a significantly high S/N to prevent adding noise when spectra are corrected for the flat field during the extraction process. For the Sun, however, as the median number of observations per day is 56, the S/N of the extracted flat field is a limitation when combining spectra reduced with the same calibrations (i.e., spectra taken on the same day). This is, for example, seen in an analysis performed in Cretignier et al. (2020, see Fig. 8 and related discussion). We therefore decided to create master flat fields that include more than five raw frames to increase the S/N and thus be able to work with daily averages of spectra without being significantly affected by the noise present in the flat field products. This implies using frames from different days, which assumes that the instrument flat field does not change significantly over short periods of time.

To test the stability of flat fields over time, we divided all the extracted flat fields by a reference one, in our case HARPN.2016-11-20T17:59:35.874_FSPECTRUM_A.fits, and analysed if the residuals were following a Gaussian statistic. To do so, we measured the ratio between the rms and the median error of the residuals for each order. We note that the errors were propagated from the FLAT calibrations. Looking at the top panel of Fig. 4, where we show this ratio, we see that this observable varies slowly over time, with sudden jumps that are due to periodic warm-ups of the detector, highlighted by vertical dashed lines (see Sect. 2.6 for more details about these warm-ups). Between warm-ups, the variation of the flat field over time is small and it is therefore possible to increase the S/N by accumulating frames over several days. However, we do not want to combine frames on both sides of a detector warm-up. Therefore, we used the following strategy to select frames for a given master flat field. To build a master flat field for date YYYY-MM-DD, we looked at all the good raw frames within a time window of ±20 days. If there is a detector warm-up closer than 20 days, for example at +10 days, we shifted the window accordingly to still have a full window of 40 days. We thus consider the raw frames between dates YYYY-MM-DD−30 days and YYYY-MM-DD+10 days. In this way, we always keep a high number of frames in master flat fields, with the drawback that master flat fields closer than 20 days from a warm-up are all identical. This is however not a problem due to the slow variation in the instrument flat field.

Fig. 4. Top: variation of the flat fields over time. Here, we show the ratio between the rms and the median error of the residuals obtained by dividing all the flats by one of reference (HARPN.2016-11-20T17:59:35.874_FSPECTRUM_A.fits). Sudden jumps can be explained by periodic detector warm-ups (dashed vertical lines). Bottom: number of raw frames used to build the master flat fields. This number has to be compared to five, which is the number of raw frames included in the old HARPS-N flat field products, and 56, which is the median number of solar observations per day.

The raw frames used for flat fielding on HARPS-N have fibres that are both illuminated by the tungsten lamp. Therefore, the flat field from the science fibre is contaminated by ghosts from the reference fibre, which are due to the secondary reflection of the echelle orders, which is not necessarily the case of a science observation with no flux, sky, thorium-argon, or Fabry-Pérot spectrum on the reference fibre. It is thus better to use calibrations where the tungsten lamp only illuminates one fibre at a time, which prevents having the ghosts from the other fibre (technique adopted on ESPRESSO). On HARPS-N, such calibrations exist and are generally used for order localisation, but nothing prevents us from using those raw frames to derive the flat field. The only problem is that there are only two such frames per calibration set (two for the science and two for the reference fibres), compared to the five frames available for standard flat fielding. However, this can be compensated for by using frames from several different days. In the lower panel of Fig. 4, we show the number of localisation raw frames included in each master flat field. As we can see, the median number of frames is 46, with the highest and lowest values at 28 and 73, which is similar to the median number of solar frames taken every day (56). Thus, when using the master flat fields described above to reduce solar data, the noise coming from flat fielding does not significantly affect daily averages of spectra. We note that all the data released with this paper have been reduced using those master flat fields.

2.5. Use of afternoon calibrations

In the old HARPS-N DRS, science data are reduced using the closest good quality calibrations back in time. This makes sense for night-time observations, as calibrations are performed in the late afternoon. However, for daily observations of the Sun, this implies that the calibrations used are from the afternoon of the previous day, therefore, at minimum more than 12 h apart. This is not critical for order localisation, flat fielding (see Sect. 2.4), background, and contamination corrections, however, it can be problematic for wavelength solutions, which set the RV zero point of the instrument every day and therefore directly impact RV stability. For each solar observation, with the Fabry-Pérot spectrum on the reference fibre, the wavelength solution is derived by taking the wavelength solution calibration (thorium-argon lamp illuminating the science fibre and Fabry-Pérot spectrum on the reference fibre) and shifting it by the Fabry-Pérot drift observed on the reference fibre. However, this assumes that the Fabry-Pérot spectrum is stable over time, which is not necessarily the case. The first HARPS-N Fabry-Pérot that was in operation from the beginning of HARPS-N in 2012–2017 November 11 showed internal drifts of a few dozen of cm s⁻¹ per day on average, with occasional sudden jumps that could be as high as 1–2 m s⁻¹. This led to its replacement in 2017 November 11 by a more stable system free of instabilities and showing reduced daily drifts. We therefore expect HARPS-N data taken after 2017 November 11 in simultaneous Fabry-Pérot mode to present less instrumental systematics.

The residual internal drift of the Fabry-Pérot implies that, in the used calibration scheme, the wavelength solution is potentially affected by a systematic error that increases with the temporal distance between the calibration and scientific exposures. While in the old DRS RVs, the data were reduced using the calibration frames of the preceding afternoon, we reduced all the solar data presented in this paper by using the closest calibrations in the future, which are the ones obtained after the solar exposures on the late afternoon of the same day. This allowed us to better mitigate the systematics induced by the internal drift of the Fabry-Pérot. We did not want to use the closest calibration in time, since morning calibrations, which are sometimes performed, would induce a jump in the RVs observed within a day, as the wavelength solution calibration used would change from the morning one to the afternoon one around midday.

2.6. Correcting for the calcium activity index contamination

The Mount-Wilson S-index, which is a measure of the chromospheric emission in the core of the Ca II H and K spectral lines, is one of the best activity proxies for solar-type stars. On the Sun, it correlates extremely well with long-term changes in total solar irradiation induced by the presence of magnetically active regions such as faculae and sunspots (e.g., Lockwood et al. 2007, Fig. 2). For other stars, this index is used to track the rotational period on the short term (e.g., Boisse et al. 2009) and magnetic cycles in the long term (e.g., Wilson 1978; Baliunas et al. 1995; Lovis et al. 2011; Dumusque et al. 2012). To be able to compare the activity of different stars, it is essential to correct the Mount-Wilson S-index from the photospheric flux contaminating the chromospheric emission in the core of the Ca II H and K spectral lines, which is performed when computing the log($R{\prime }_{\mathrm{HK}}$) activity index (Noyes et al. 1984). This index is often used to model stellar activity in RV time series (e.g., Dumusque et al. 2011a; Díaz et al. 2016) or to train non-parametric models such as Gaussian processes to mitigate stellar activity (e.g., Rajpaul et al. 2015; Langellier et al. 2020).

However, the S-index or log($R{\prime }_{\mathrm{HK}}$) is difficult to measure since the related chromospheric emission appears in the extreme blue. At those wavelengths, more absorption from the Earth’s atmosphere and lower sensitivity of optical fibres, optics, and the detector make it difficult to reach high S/N. In addition, the chromospheric emission happens in the core of the extremely saturated Ca II H and K spectral lines, where the S/N is even lower. Therefore, a careful background correction has to be performed to extract the chromospheric emission correctly. In addition, on HARPS-N, a small leak in the detector cryostat causes humidity to build up on the detector, which increases its reflectivity over time. A periodic warm-up of the detector solves the issue for a little while. However, with the reflectivity increase, the brightness of ghosts increases as well, modifying the contamination in the extreme blue with time. On the first and second panels of Fig. 5, we show the bluest part of two raw solar science frames of similar S/N obtained with HARPS-N (frames HARPN.2017-04-09T10-06-46.964.fits and HARPN.2017-04-12T09-09-40.30.fits), with the Sun on the science fibre and the Fabry-Pérot on the reference fibre. The images were taken just before and after a warm-up of the detector, and they show how the ghosts significantly reduce in intensity with the intervention.

Fig. 5. Raw images from HARPS-N taken before (HARPN.2017-04-09T10-06-46.946, left panel) and after (HARPN.2017-04-12T09-09-40.305, middle panel) a warm-up of the detector, and their ratio (right panel). The location of the core of the Ca II K and H lines in orders one, two, and three is highlighted by grey boxes. As we can see, the core of the Ca II K line in order one (K, bottom left) and the H line in order three (H) are strongly contaminated by ghosts. In addition, we also see that the contamination due to ghosts is significantly reduced with the warm-up process, which removes the humidity accumulated over time on the detector due to a small leak in its cryostat.

In Fig. 5, we highlighted the location of the Ca II H and K lines. We note that the core of the Ca II K line appears on two consecutive orders, extracted orders one and two, and the core of the Ca II H lines appears on the third extracted order. The core of the K line in order one and the H line in order three are significantly affected by ghosts. With humidity building up on the detector over time, ghost contamination increases, and therefore it increases the measured log($R{\prime }_{\mathrm{HK}}$). After a warm-up, the contamination is significantly reduced, which induces an abrupt offset in the log($R{\prime }_{\mathrm{HK}}$) time series. This is clearly visible in the top left panel of Fig. 6, where we show the raw log($R{\prime }_{\mathrm{HK}}$) extracted by the old HARPS-N DRS along with the time localisation of each warm-up.

Fig. 6. Top left: comparison between the raw (dots with yellow to blue colour scale) and corrected (red dots) log($R{\prime }_{\mathrm{HK}}$) calcium activity index, as described in the text. Bottom left: differences between the raw and corrected log($R{\prime }_{\mathrm{HK}}$) time series. We clearly see that the slow increase in the activity index followed by sudden jumps at each warm-up of the detector seen in the raw log($R{\prime }_{\mathrm{HK}}$) time series is strongly reduced after correction. Top and bottom right: closer look at the regions around the third and fourth warm-ups.

The old and new HARPS-N DRS do not correct for ghost contamination, as this is challenging to perform. First of all, ghosts are at low S/N which makes their extraction difficult. In addition, ghosts are oblique relative to echelle orders and they contaminate them over hundreds of pixels. Finally ghosts, which are due to the secondary reflection of the echelle orders, are lower flux copies of the solar spectrum and thus do not present a flux that is smoothly varying along the direction of the ghost. Because it is difficult to know which ghost corresponds to which echelle order, we decided to use a rather simple technique, described in Appendix D, to correct for the contamination.

The log($R{\prime }_{\mathrm{HK}}$) measured with and without the correction can be seen in the top left panel of Fig. 6 and the difference between the two is shown in the bottom left panel. It is clear that the increase in flux over time followed by sudden jumps at each detector warm-up is strongly mitigated in the corrected time series. A closer look at the regions around the third and fourth warm-ups can be seen in the right panel of the same figure. For the third warm-up, the observed offset without correction is 0.015 dex, and it goes down to 0.003 dex after correction. For the fourth warm-up, the offset goes from 0.024 to 0.002 dex. The correction applied therefore strongly mitigates the effect of ghost contamination in the log($R{\prime }_{\mathrm{HK}}$) time series, and it is difficult to say at this stage if the tiny offset remaining at the third and fourth warm-ups are from solar or from instrumental origin.

The corrected log($R{\prime }_{\mathrm{HK}}$) activity index that we can see in the top left panel of Fig. 6 has an exquisite precision due to the high S/N solar observations. The median error bar is 0.0014 dex for 5-min integration times, with a standard deviation of 0.0002 dex. As we can see, the 25-day rotational modulation due to active regions turning with the Sun is clearly visible. Detailed analyses of this log($R{\prime }_{\mathrm{HK}}$) time series were performed in Maldonado et al. (2019) and Milbourne et al. (2019). We note that in those papers, the log($R{\prime }_{\mathrm{HK}}$) time series were corrected with a similar method as described here.

3. The HARPS-N solar data set

3.1. Data selection

For this paper, and the related data release, we only selected good quality spectra over the first three years of operation of the HARPS-N solar telescope. As the Sun is resolved with the 3-inch entrance lens of the solar telescope, any kind of absorption, such as clouds, not evenly distributed over the solar disc induces a huge Rossiter–McLaughlin-type effect, with RVs that can depart from the true value by several hundreds of meters-per-second, as the v sin i of the Sun is 2 km s⁻¹. To determine if spectra are contaminated by clouds, Collier Cameron et al. (2019) performed a daily regression of −2.5log[(S/N)²] in echelle order 60 (as a proxy for apparent solar magnitude) against airmass. Outliers from the expected linear extinction relation each day were considered to have been affected by clouds, and they were assigned probabilities of belonging to a separate outlier population using a Gaussian mixture model (Foreman-Mackey et al. 2014; Hogg et al. 2010). For the present paper, and the related data release, we only considered the data for which the probability of being of a good quality is ≥0.99.

With such a threshold, 65% (34 550) of the spectra taken with the HARPS-N solar feed are selected out of the 53 417 available. The median cadence between consecutive observations, excluding the night gaps, is 5.42 min, which is consistent with the effective time between measurements, being 5 min exposure time plus 20 s readout. We note that we chose an exposure of 5 min to average out the signals induced by p-mode oscillations (e.g., Chaplin et al. 2019).

3.2. Comparison between the solar radial velocities reduced with the new and old HARPS-N DRS

In this section, we compare the three years of solar RVs already published in Collier Cameron et al. (2019), which were obtained with the old HARPS-N DRS and referred to as the old solar RVs, with the RVs obtained with the ESPRESSO DRS and the improvements discussed in this paper (i.e., the new HARPS-N DRS), referred to as the new solar RVs. When estimating RV errors, the old HARPS-N DRS included photon noise, drift noise induced by the drift measurement performed on the simultaneous Fabry-Pérot spectrum, and the error in wavelength calibration, extracted from the polynomial fit performed to derive it. The median value for this error in wavelength calibration is 0.35 m s⁻¹, which is much less than the 0.74 m s⁻¹ measured in Sect. 2.2. Thus, the HARPS-N RVs derived using the old DRS have error bars that are significantly underestimated.

In the new DRS, the estimated RV errors only include photon and drift noises, because to properly account for the wavelength solution noise, which impact all the observations performed within the same calibration set in the same way (obtained on average every day), we need to consider a covariance matrix which is block-diagonal; each block regroups all the observations attached to the same wavelength solution (Delisle et al. 2020). Each block should thus have the same calibration noise of 0.49 m s⁻¹ for the new DRS data. Another way to include this calibration noise correctly is to daily-bin the RV data and then add a calibration noise of 0.49 m s⁻¹ to the RV errors in quadrature.

In the second panel of Fig. 7, we show the comparison between the old and new solar RVs. The colour scale for the new RVs corresponds to the S/N of the measurements. As we can see, the S/N degraded from ∼460 to ∼350 in three years, thus a reduction in flux of nearly a factor of two. We never changed the exposure time nor the neutral density filter used to attenuate the incoming flux, thus this could be due to ageing of the solar telescope plexiglass dome, the optics of the solar telescope, or the fibre guiding solar light into HARPS-N or a mix of the three. The neutral density filter was changed recently to solve this issue. We note that around BJD = 2458100 (December 2017), the fibre guiding the light from the solar telescope to HARPS-N was damaged following an intervention on the HARPS-N calibration unit. This led to observations with very low S/N, giving RVs that were not precise enough for this study. Those data were rejected from the data released, explaining the absence of data around BJD = 2458100.

Fig. 7. Comparison between the solar RVs reduced using the old HARPS-N (blue) and the new HARPS-N DRS (orange). Top panel: weighted rms for each chunk of RV data of 50 days, which corresponds to two solar rotations and therefore include activity signal. Orange and blue bars correspond to the weighted rms measured on the new RVs and on the old RVs, respectively. The number for each chunk corresponds to the improvement in weighted rms brought by new DRS (in percentage). Second panel: solar RVs as derived by the old and new HARPS-N DRS. The blue and orange lines correspond to a smooth version of the calcium activity index variation linearly fitted to the RVs, to remove the strong long-term correlation between RVs and the log($R{\prime }_{\mathrm{HK}}$), and thus to mitigate the contribution from the solar magnetic cycle (see text). The colour scale for the new RVs corresponds to the S/N of the measurements. We note that we lost a significant amount of data around BJD = 2458100 (December 2017) due to a damaged fibre (see text). Third panel: difference between the new and old solar RVs, for which a long-term trend is clearly visible. We highlighted with red and green vertical lines the significant discrepancies in those RV differences. The red lines correspond to problems in the old DRS RVs, while for the green vertical line, it is not clear if the correct RVs are the ones derived with the old or new DRS. Bottom panel: generalised Lomb-Scargle periodogram of the residuals after removing the long-term trend observed in the RVs. The small inset shows a zoom around the periods 10–30 days to highlight stellar rotation.

To see the improvements brought by the new DRS, we divided the data set in chunks of 50 days, which corresponds to approximately two solar rotations, and measured the weighted rms within each chunk (values reported in the top panel of Fig. 7). We find that 17 of the 21 chunks show a greatly reduced weighted rms in the new solar RVs, with an increased precision of up to 24%. This improvement was mainly brought on by (i) the new drift map method with tailored thorium line selection to derive wavelength solutions (see Sect. 2.2), which reduces calibration to calibration RV offsets by 0.55 m s⁻¹, and by (ii) the use of the afternoon calibrations rather than the ones from the previous day. In the left panel of Fig. 8, we show the day-to-day RV offsets measured on the old and new RVs after daily-binning the data. As we can see, the new RVs present a distribution with a much smaller MAD of 1.07 m s⁻¹, compared to 1.27 m s⁻¹ for the old RVs. Thus, when analysing day-to-day RV offsets, the new solar RVs show a very significant improvement of $\sqrt{1.{27}^2-1.{07}^2}=0.68$ m s⁻¹.

Fig. 8. Left: histogram of the day-to-day RV offsets measured on the old and new solar RVs after daily-binning the data. Right: intra-day RV variation obtained by removing the daily mean on each single day.

In the right panel of Fig. 8, we show the RV variation on timescales shorter than a day, which is obtained by subtracting the daily mean of each day. On this short timescale, we measured a very similar MAD for the new and old DRS, 0.54 and 0.55 m s⁻¹, respectively, which is not surprising as all the data within a day generally use the same set of calibrations. This MAD is significantly larger than the median photon noise of 0.25 m s⁻¹. This discrepancy is likely due to stellar granulation (e.g., Meunier et al. 2017a; Dumusque et al. 2011b) and oscillations that are not fully averaged out with an integration time of 5 min (e.g., Chaplin et al. 2019), as we know that those stellar signals are of several dozens of meters-per-second on the Sun (e.g., Lefebvre et al. 2008).

In the third panel of Fig. 7, we can see the difference in RVs between the new and old solar RVs. We see that the two data sets follow each other closely, except for a few points highlighted with red and green vertical lines. After visually inspecting in the old and new DRS RVs if the corresponding data points were following the general RV variation or if they were outliers, we concluded that the 12 red vertical lines correspond to discrepancies induced by the old DRS. The one green line corresponds to data for which it is difficult to say if the old or the new DRS gives a better solution.

Another significant aspect that can be seen in the RV differences is a small long-term drift that is rather linear with time. Except probably for the behaviour observed after BJD = 2458070, this drift does not correlate with the difference in wavelength solutions between the old and new DRS well (see right panel of Fig. 3). Therefore, it does not seem that it comes from the way wavelength solutions are derived. One thing that is different between the old and new DRS is the computation of the cross-correlation function (CCF, Baranne et al. 1996). First, the old reduction uses the HARPS-N G2 template to perform the cross-correlation, while the new DRS uses the ESPRESSO G2 template. Those templates consist of the wavelength of all the spectral lines used to compute the CCF, with some weights corresponding to the strength of the lines. As described in Pepe et al. (2002), lines that are stronger have a better RV precision and therefore should be weighted more when computing the CCF. In the old DRS, the CCF is computed using the normalised depth of the lines with respect to the continuum as weight, as described in Pepe et al. (2002), while the new DRS uses the square of the line depth, which corresponds to the true RV precision of a spectral line better (Lovis et al. 2010). As lines in the blue part of the spectrum of a G2 star are very strong and their density is higher than in the red, using the depth of the lines squared as weight put more importance on the blue part of the spectrum, as shown in Fig. 9. Several studies demonstrate on simulations that the blue part of the visible is more sensitive to stellar activity (e.g., Reiners et al. 2010; Desort et al. 2007). This has also been shown on real data, but generally on stars more active than the Sun (Zechmeister et al. 2018; Feng et al. 2017; Anglada-Escudé & Butler 2012). If it is the case for the Sun, we expect the new solar RVs to be more sensitive to stellar activity than the old ones.

Fig. 9. Comparison of the weights for all the lines used to build the CCF between the old and new DRS G2 templates.

On stars such as the Sun, the main source of activity is induced by faculae (e.g., Meunier et al. 2010a; Haywood et al. 2016; Milbourne et al. 2019). The magnetic field in faculae inhibits locally stellar convection, which suppress the convective blueshift locally and therefore induces a net redshift (e.g., Lindegren & Dravins 2003; Meunier et al. 2010a). Thus, it is expected that the calcium activity index, probing faculae, and the RV exhibit a strong correlation (e.g., Lindegren & Dravins 2003; Lovis et al. 2011; Dumusque et al. 2011a; Meunier et al. 2017b). By looking at the correlation between the RV and the log($R{\prime }_{\mathrm{HK}}$) time series in Fig. 10, we clearly see that in the long term, both the old and new solar RVs show a strong correlation as expected. Furthermore, the new solar RVs present a higher Pearson correlation coefficient of R = 0.93 compared to R = 0.77 for the old RVs. Regarding short-term variations, the new solar RVs correlate slightly better, with R = 0.57 compared to R = 0.51, however, the difference is less significant. This analysis confirms that the new solar RVs are more sensitive to stellar activity than the old ones, and as explained in the preceding paragraph, this is likely because the new DRS uses a CCF template that puts more weight on the blue part of the spectrum.

Fig. 10. Top: correlation between the RVs and the log($R{\prime }_{\mathrm{HK}}$) for long-term variations, which in this case means that the time series were smoothed using a 30-day rolling average. Bottom: correlation between the short-term variations, meaning the raw time series minus the smoothed ones. In the legend, we show the values for the R Pearson correlation coefficient for the new and old solar RVs.

In the bottom panel of Fig. 7, we compared the Generalised Lomb-Scargle periodograms (Zechmeister & Kürster 2009) of the new and old solar RVs. As it seems that the long-term drift is induced by the solar magnetic cycle, we removed its contribution by linearly fitting a smooth version of the calcium activity index variation⁷ to the RVs, as performed in several studies (e.g., Delisle et al. 2018; Díaz et al. 2016). The best-fitted linear models are represented as smooth lines in the second panel of the same figure. At first glance, the two periodograms are nearly identical. As we can see, the most significant signals can be found around 27 and 13.5 days, the synodic rotational period of the Sun, and its first harmonic, which is expected from stellar activity (e.g., Saar & Donahue 1997; Boisse et al. 2009). In the small inset in the bottom panel of the figure, we clearly see that those signals modulated by stellar rotation have more power in the new DRS RVs, mainly for the peak at 13.4 days. This can be explained by the fact that the new solar RVs are less affected by instrumental and data reduction systematics and, therefore, the signal from stellar activity is better characterised. We note, however, a slight excess of power around 180 days in the new DRS RVs. This periodicity is also observed in independent observations of the Sun over magnetic cycle 24, estimated using images from the Helioseismic and Magnetic Imager onboard the Solar Dynamics Observatory, in the RV variations as well as in the plage filling factor and unsigned magnetic flux (see Fig. 5 in Haywood et al. 2020). A periodicity of ∼155 days is also found in the sunspot blocking function, the 10.7 cm radio flux, and the sunspot number of solar cycle 19, 20, and 21 (Lean & Brueckner 1989) and in the recurrence of hard solar flares in cycle 21 (Rieger et al. 1984). We conclude that this 180-day periodicity is probably of a stellar origin and not from instrumental systematics, however further investigation is required.

3.3. Data release

Alongside this paper, we release the first three years of HARPS-N solar data. Those data consist of the following data products reduced with the new HARPS-N pipeline:

– The extracted echelle-order spectra, corrected from the instrumental blaze, in the Earth rest-frame. These products are called S2D spectra due to their two-dimensional shape.

– The extracted merged spectra, corrected from the instrumental blaze, in the Solar System barycentric rest-frame. These products are called S1D spectra.

– The CCFs obtained by cross-correlating the S2D spectra with a synthetic mask derived from solar observations (ESPRESSO G2 mask), thus optimised for the Sun.

In addition to those data products, we also release the time series of relevant quantities such as the RV, the BIS SPAN, the FWHM, and different activity indicators. They are all listed in Table A.1 with their descriptions.

The different data products as well as the time series can be accessed using the Data Analysis Center for Exoplanet (DACE) web platform⁸, or using the DACE python API⁹ (with the solar spectroscopy module). More detailed information about the available data and how to access them can be found in Appendix A.

4. Discussion and conclusion

In this paper, we describe a new data reduction for HARPS-N that allowed us to obtain unprecedented RV precision for the three years of solar data that are released alongside this manuscript (see Sect. 3.3 and Appendix A for more information about the released products). We have also compared the newly derived RVs with the ones already published in Collier Cameron et al. (2019) and clearly demonstrated that the new data are less affected by instrumental and data reduction systematics.

Compared to the Collier Cameron et al. (2019) solar RVs data that were reduced using the current HARPS-N DRS (version 3.7, referred to as the old DRS), we reduced new RVs using the ESPRESSO DRS tailored to HARPS-N (version 2.2.3, referred to as the new DRS, see Sect. 2.1). In addition, we developed a new method to derive wavelength solutions (see Sect. 2.2 and Appendix B) and curated the list of thorium lines that are used when deriving a wavelength solution, so that only the most stable thorium lines are used over time (see Sect. 2.3 and Appendix C). We also reduced the data taking the closest calibrations in the future rather than in the past (see Sect. 2.5) to mitigate any systematics induced by Fabry-Pérot internal drifts, and we used master flat fields that have a much higher S/N than standard flat fields to prevent the solar observations being limited by them (see Sect. 2.4). We also extracted a calcium activity index, which was corrected from ghost contamination (see Sect. 2.6 and Appendix D).

We performed several analyses in Sect. 3.2 to compare the new and old DRS solar RVs. When dividing the solar data in chunks of 50 days to still see the 25-day effect of stellar rotation, we observe that for most of the chunks, the new solar RVs are less affected by noise, with an improvement of up to 24% (see Fig. 7). At first, we could think that the new DRS mitigates part of the stellar activity signal. We see below when discussing the correlation between the RVs and the activity index log(R${\prime }_{\mathit{HK}}$) that this does not seem to be the case. In fact, a similar improvement is seen when studying the RV variability on a timescale of a day, which is more affected by instrumental and data reduction systematics as calibrations are performed on a daily basis. On that timescale, it is clear from the left panel of Fig. 8 that the new RVs are less affected by noise since the measured MAD goes down to 1.07 m s⁻¹ with the new DRS, compared to 1.27 m s⁻¹ with the old one. As systematics adds up in quadrature, this implies that the new reduction reduces a noise of 0.68 m s⁻¹ present in the old solar RVs. This improvement is mainly brought by the new wavelength solution derivation method, the optimal selection of thorium lines, and the use of calibrations closer in time.

On a daily basis, a noise floor of 1.07 m s⁻¹ is still much larger than the photon noise error of a few centimetres-per-second¹⁰ and the wavelength calibration noise error, estimated to be 0.49 m s⁻¹ when analysing the RV offsets between consecutive wavelength solutions (see Sect. 2.2 and Fig. 2). Although we cannot exclude some remaining instrumental and data reduction systematics, it is likely that the remaining 0.96 m s⁻¹ ($\sqrt{1.{07}^2-0.{03}^2-0.{49}^2)}$) in RV variation over a day are induced by stellar signals. On that timescale, several studies have shown that super granulation on solar-type stars can induce an RV variability at the level of a meter-per-second (e.g., Lefebvre et al. 2008; Dumusque et al. 2011b; Meunier et al. 2015).

We also looked at the RV variability within a day. In the right panel of Fig. 8, we see that the old and new data reductions give a similar MAD of 0.55 and 0.54 m s⁻¹, respectively. As all data within a day are generally reduced with a single calibration set, and assuming that the Fabry-Pérot étalon that is used to correct for instrumental drift does not drift internally, we do expect to measure an RV variability close to photon noise. As for a single exposure, the median photon noise is 0.25 m s⁻¹, and the solar data are clearly affected by some extra variability. This noise is likely due to stellar granulation (e.g., Meunier et al. 2017a; Dumusque et al. 2011b) and oscillations that are not fully averaged out with an integration time of 5 min (e.g., Chaplin et al. 2019), as we know that those stellar signals are of several dozens of meters-per-second on the Sun (e.g., Lefebvre et al. 2008).

We observed that the new solar RVs show a stronger drift with time than the old RVs (see Fig. 7), and this drift cannot be explained by the use of the new selection of thorium lines and drift map wavelength solutions (see Sect. 3.2). The presence of this stronger drift induces a stronger correlation of R = 0.93, compared to R = 0.77, between the RVs and the log($R{\prime }_{\mathrm{HK}}$) in the long term (see Fig. 10), which is expected in the case of a magnetic cycle as seen on the Sun (e.g., Dumusque et al. 2011a; Meunier et al. 2010b; Lindegren & Dravins 2003). After investigation, it seems that this drift is mainly induced by the use in the new DRS of a new cross-correlation template, the ESPRESSO G2 mask, that puts more weight on the blue part of the spectrum. The blue is more sensitive to stellar activity effects (e.g., Zechmeister et al. 2018; Feng et al. 2017; Anglada-Escudé & Butler 2012; Reiners et al. 2010; Desort et al. 2007), and thus the RVs are more affected by the solar magnetic cycle. We believe that being more sensitive to stellar activity can be seen as an improvement, as this increases the S/N of the activity signal, which helps techniques that needs to better characterise stellar activity before mitigating it (e.g., Rajpaul et al. 2015; Cretignier et al. 2020; Collier Cameron et al. 2020; de Beurs et al. 2020). However, we have to keep in mind that this complicates the search for planetary signals if no stellar activity correction is performed.

On the timescale of solar rotation, the old and new DRS gives a similar correlation of R = 0.51 and R = 0.57. This confirms that although the log($R{\prime }_{\mathrm{HK}}$) can be used to efficiently model the RV variations induced by solar-like magnetic cycles (e.g., Lovis et al. 2011; Meunier & Lagrange 2020), it cannot be used to model the RV stellar activity signals on the rotational timescale efficiently.

From the different analysis performed in this paper, it is clear that the new solar RVs and corresponding high-resolution spectra and CCFs should be preferred compared to the old DRS products, as the instrumental and data reduction RV variability is smaller, and also because those data better reflect the effect from the solar magnetic cycle. After releasing these data through the DACE web platform (see end of Sect. 3), the next effort from the HARPS-N collaboration will be to run the new DRS with all the described improvements in online mode at the Telescopio Nazionale Galileo and to re-reduce all the HARPS-N data, as the gain in terms of precision and stability is clearly demonstrated on the Sun.

On a final note, at best, we reduced all the observed systematics on the HARPS-N solar data set. The release of this unprecedented time series of spectra and radial velocities is crucial to improve our understanding of stellar signals affecting the Sun and solar-type stars. We hope that the community will use such data to develop novel methods to mitigate stellar signals in radial-velocity data sets, with the goal of enabling the detection of other Earths.

Acknowledgments

XD is grateful to The Branco Weiss Fellowship-Society in Science for its financial support. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 851555/SCORE). ACC acknowledges support from the Science and Technology Facilities Council (STFC) consolidated grant number ST/R000824/1. This work has been carried out in the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). AM acknowledges support from the senior Kavli Institute Fellowships. This research used the DACE platform developed in the frame of PlanetS (https://dace.unige.ch). This work was performed under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute (RDH).

References

Appendix A: Description of the data released

The data release consists of the following data products:

– The extracted echelle-order spectra, corrected from the instrumental blaze, in the Earth rest-frame. These products are called S2D spectra due to their two-dimensional shape. The first and second extensions in the FITS file contain the blaze-corrected extracted flux per pixel and corresponding error for each spectral order. The error corresponds to the photon plus read-out noises of the detector added in quadrature, and divided by the blaze, so that the corresponding error can be directly used with the flux given in the first extension of the FITS file. The third extension corresponds to the quality flag of the pixels for each order, zero being good and anything else being bad. Hot and bad pixels are flagged that way. Extensions four and five are the wavelength solution in the vacuum and in the air, and extensions six and seven are the width of pixels in wavelength in the vacuum and in the air, respectively. We note that all wavelengths are in Ångström. Because of dispersion, the size of each pixel in Ångström changes with wavelength, which implies that for a given order, the continuum of an S2D spectrum shows a significant slope. To get a flat continuum, the easiest way is to transform the wavelengths in a logarithmic scale.

– The extracted merged spectra, corrected from the instrumental blaze, in the Solar System barycentric rest-frame. These products are called S1D spectra. The only extension in the FITS file includes the wavelength in the vacuum and in the air for each point of the merged spectrum, its flux, and the quality flag of the point, as defined in the first item above. We note that merged-1D spectra are interpolated on a grid constant in velocity space and not in wavelength space. The step between each point is 0.82 km s⁻¹, equivalent to the width of a pixel in velocity space.

– The CCFs obtained by cross-correlating the S2D spectra with a synthetic mask derived from solar observations, thus optimised for the Sun. The first extension in the FITS file gives the CCF measured for each echelle order, with a step of 0.82 km s⁻¹, in addition to the photon noise weighted average CCF over all orders. Therefore, the extension has the shape N_CCF × (N_ord + 1), where N_CCF is the number of points of the CCF and N_ord is the number of echelle orders, 69 for HARPS-N. The second extension gives the photon noise plus read-out errors, and the third one gives the quality flag of each point as defined in the first item above.

In addition to those data products, we also release the time series of relevant quantities such as the RV, the BIS SPAN, the FWHM, and different activity indicators. They are all listed in Table A.1 with their descriptions.

The data can be accessed using the Data Analysis Center for Exoplanet (DACE) web platform¹¹, or using the DACE python API¹² (with the solar spectroscopy module). Bulk downloads are available on the main web page to download all the available products (S2D, S1D, and CCF), only the CCF products, or the time series. Using the solar spectroscopy data base (go to the main web page and then follow the link to the Solar spectroscopy database), it is possible to scan through all the released data, select specific observations by using filters (magnifying glass on the right of each column) and download the S2D, S1D, or CCF products of specific observations by selecting the rows of interest (click on them), and then using one of the S2D, S1D, or CCF buttons in the top right corner of the table. We note that although by default only eight columns are shown on the Solar spectroscopy database web page, it is possible to get access to more information by clicking on the gearwheel in the top left corner of the table. All this can also be done using the DACE python API for people who want to directly access the data within a python script¹³. A small tutorial to perform a basic analysis on the solar data released is also available online¹⁴.

Appendix B: Wavelength solution derived using the drift map method

In Fig. B.1 we show the different steps performed to derive and validate the quality of a wavelength solution using the drift map method. The first step is to measure the drift of the thorium lines between a thorium calibration, for which we want a wavelength solution, and a reference thorium calibration, for which a reference wavelength solution already exists. For the work presented in this paper, we derived a reference wavelength solution for 2017 November 11 by combining information from a thorium calibration followed by several consecutive Fabry-Pérot calibrations, which allowed us to obtain a very high-S/N Fabry-Pérot spectrum and, therefore, a precise reference wavelength solution. On the measured drift of all thorium lines, we then fitted a smooth two-dimensional polynomial of shape ΔX = ∑_i, jXⁱO^j, where i = 0, 1, 2 and j = 0, 1, 2, 3, 4, X corresponds to the pixel position in the dispersion direction, and O to order number, to account for variations in the cross-dispersion direction. The corresponding wavelength solution was obtained by shifting the reference wavelength solution using the fitted 15 free parameters’ two-dimensional polynomial. We note that we selected a two-dimensional polynomial with i ≤ 2 and j ≤ 4, as tests using polynomials of smaller or higher dimensions gave residuals with a larger scatter or no significant improvement, respectively.

Fig. B.1. Derivation of the wavelength solution for the thorium calibration HARPN.2017-03-25T16-47-52.890.fits using the drift map method with the frame HARPN.2017-11-12T15:39:58.405.fits selected as a reference. Top: first, we measured the drift of all the thorium lines between those two calibrations. Each thorium line is represented as a dot with a size proportional to the error in drift and colour to the drift measured. We note that the thorium lines have been represented at their exact location on the detector using the localisation of the orders (grey dotted lines). Bottom left: then, we fitted a smooth two-dimensional polynomial on the measured drift of all the thorium lines. The corresponding wavelength solution was obtained by shifting the reference wavelength solution with the fitted two-dimensional polynomial. Bottom right: measured drift on the thorium lines between the thorium calibration and the reference before adjusting the two-dimensional polynomial (top) and after fitting for it (bottom). As we see, the χ2 of the two-dimensional polynomial fit is close to one, and the residuals show a scatter, which was measured with the MAD, which is significantly reduced.

Appendix

Tailored selection of thorium lines for the new HARPS-N DRS

With the drift map method used to derive wavelength solutions, described in Sect. 2.2 and Appendix B, we only used 15 free parameters. Thus, we can significantly reduce the number of thorium lines to be used, as long as we are left with thorium lines distributed all over the detector. In addition, the new DRS only fits for a single Gaussian component for each main thorium line, while the old DRS fitted for several components when blends were known (Lovis & Pepe 2007). If blends vary with time in a different way than the main component, we observe an asymmetric variation in the main component, which modifies the line’s barycentre and therefore induces a spurious RV drift. This spurious change affects the derived wavelength solution and, therefore, the RV stability of the spectrograph in the long term. As a result, we studied the behaviour of all the thorium lines used in the drift map method and rejected the ones that showed strange behaviours.

As described in Sect. 2.2, due to ageing of the thorium-argon lamp used for wavelength solution, the number of lines used over time decreases quite significantly. Starting from the line selection of Lovis & Pepe (2007), we selected the lines that were used 99% of the time over the all HARPS-N lifetime¹⁵. At the end of this process, we were left with 3061 thorium lines.

Then we looked at how those remaining lines were behaving as a function of time compared to the median of all the lines. We observed that some lines showed some drifts over time of dozens of meters-per-second-per-year compared to the median of all the lines, and some others were strongly correlated with the thorium flux ratio, which is the total flux emitted by the thorium lamp at a given time relative to a reference time. The flux ratio is a good proxy to track the ageing of the thorium lamp over time. In the top panel of Fig. C.1, we show the thorium flux ratio of the thorium lamp used on HARPS-N from the commissioning in 2012–June 2020. As we can see, the thorium flux ratio changes significantly over time, mainly after BJD = 2458070, where we observe a strong linear drift, characteristic of a thorium lamp ageing fast. In the middle and bottom panel of the same figure, we show the RV drift of a thorium line that shows a strong correlation of R = 0.80 with the thorium flux ratio (in the old DRS). This is likely explained by this thorium line being affected by a blend which varies with the ageing of the lamp. It is clear that such a thorium line shows a strong systematics over time and should not be used.

Fig. C.1. Top: thorium flux ratio as a function of time, which is the total flux emitted by the thorium lamp at a given time relative to a reference time. Middle: RV drift of the thorium line at 4833.77 Å measured with the old (blue) and new (orange) DRS. On both data sets, we fitted a model corresponding to a linear drift plus a linear correlation with the thorium flux ratio. In the legend, we report the measured slope of the linear drift. Bottom: RV drift residuals after subtracting the fitted linear drift. In the legend, we report the R Pearson correlation coefficient between those residuals and the thorium flux ratio, as well as the ratio between the rms of those residuals and the median error in RV for this spectral line to detect a non-Gaussian behaviour.

To assess the stability of the remaining 3061 thorium lines that are almost always used, we measured three different quantities. For each line, we fitted a model corresponding to a linear drift plus a linear correlation with the thorium flux ratio. We then considered the slope of the linear drift and the Pearson correlation coefficient between the RV drift residuals after removing the fitted slope and the thorium flux ratio as two relevant quantities of line stability. We also measured the ratio between the rms of the RV drift residuals and the median RV error, which corresponds to our third quantity. This last quantity probes the presence of remaining signals in the RV drift residuals. We measured those quantities for all the thorium lines used in the old and new DRS, and we report their values in Fig. C.2.

Fig. C.2. Left: slope measured on the RV drift of all the thorium lines used in the old (blue) and new (orange) DRS. Middle: Pearson correlation coefficient between the RV drift of each thorium line and the thorium flux ratio (see Fig. C.1). Right: rms of the RV drift residuals when subtracting the measured linear drift, divided by the median error in RV. For each subplot, the red vertical dashed lines correspond to the cut-offs used to reject thorium lines showing strong systematics (±10 m s−1 yr−1 for slope, ±0.3 for Pearson correlation coefficient, and < 1.5 for RV residual rms over error in RV). The black histogram corresponds to the selection of lines used 99% of the time, and the orange one corresponds to the final selection used in the new DRS after rejecting lines using the different cut-offs.

Looking at the different distributions obtained in Fig. C.2 for the new DRS thorium line selection, we decided to reject the lines for which the slope was larger than ±10 m s⁻¹ yr⁻¹, for which the Pearson correlation coefficient with the thorium flux ratio was larger than 0.3, and for which the residual rms over the median error in RV was larger than 1.5. We chose those cut-offs to reject lines with the strongest systematics, while keeping a minimum of nine thorium lines per echelle order to sufficiently constrain the two-dimension polynomial fit performed to derive wavelength solutions (see Appendix B). After this drastic selection, we were left with 2081 thorium lines.

Finally, as we used the Lovis & Pepe (2007) catalogue as a starting point, the thorium line wavelengths come from Palmer & Engleman (1983). Redman et al. (2014) presented a new atlas of thorium lines with a more precise wavelength determination and uncertainties. We therefore cross-matched our final selection of lines with this updated catalogue to provide the most up-to-date wavelengths and uncertainties for thorium lines. To do so, we required that a line in our selection had strictly one line in the Redman et al. (2014) atlas within 1 km s⁻¹. We rejected the 37 lines that did not match this criterion. Finally, we checked that the ionisation level between the two catalogues was the same. This was not the case for three lines, and therefore we rejected them as well. Our final selection of thorium lines, with updated Redman et al. (2014) wavelengths and uncertainties, is thus composed of 2041 lines that can be found in Table C.1. We can also see the distribution of this final selection of lines over the HARPS-N detector in Fig. C.3.

Fig. C.3. Position of the thorium lines used by the old (blue) and new (orange) DRS as a function of their position on the detector. We note that the position in dispersion direction is in pixels, and the position in the cross-dispersion direction is in echelle order number.

Appendix D: Correcting the calcium activity index from ghost contamination

To calculate the contamination induced by ghosts on the measured flux in the core of the Ca II H and K lines, we optimally extracted the flux, using the profile of the corresponding echelle orders, on both sides of the orders containing the Ca II H and K lines. To determine on which pixel to extract the flux, we used the order localisation and shifted the localisation by plus or minus eight pixels in the cross-dispersion direction. We highlighted the localisation of those pixels as blue and orange lines in the left panel of Fig. D.1. As the separation between the science and calibration fibre is at least sixteen pixels, these eight pixel shifts allowed us to extract the flux from ghosts, without being significantly contaminated by the solar or Fabry-Pérot spectra. Optimally extracting the flux using the profile of the orders allowed us to measure the exact flux that contaminates the science fibre. In the caption of Fig. D.1, we describe how the ghost contamination is measured in the core of the Ca II H line in detail. We performed the same measurements for the core of the Ca II K line appearing in order one and two of HARPS-N, and we finally corrected from this contamination before measuring the corrected S-index and relative log($R{\prime }_{\mathrm{HK}}$).

Fig. D.1. Derivation of the chromospheric emission in the core of the Ca II H line with estimation of the contamination induced by the perturbing ghost. Left: raw image centred on extracted order three with the position of the Ca II H line core. Right: extracted spectra at different positions. The green spectrum is the stellar spectrum extracted from order three (green line on the left panel). The blue and orange spectra are the flux optimally extracted, using the same order profile as order three, minus and plus eight pixels in the cross-dispersion direction on each side of order three (blue and orange lines on the left panel). We note that the extracted fluxes for those spectra have been multiplied by 20 to highlight the contamination from the ghost (red dashed line in the left panel). The vertical green dashed line corresponds to the centre of the Ca II H line core (green dot on the left panel), and the blue and orange vertical dashed lines correspond to the centre of the ghost on both sides of order three (blue and orange dots on the left panel). Assuming that the ghost does not show any curvature, which is a good approximation from one side of the order to the other, we measured the crossing of the ghost and order three (red dot in left panel, and red vertical dashed line in the right panel) as the middle between the orange and blue dots. To measure the chromospheric emission in the core of the Ca II H, the extracted flux from the green spectrum was multiplied with the triangular green response of width 1.09 Å (Vaughan et al. 1978). To measure the effective ghost contamination, we multiplied the contaminating flux on the blue and the orange spectra with the same triangular response, taking into account that the ghost does not cross order three at the centre of the Ca II H line core (the difference is highlighted by the red arrows). Finally, the contamination inside the core of the Ca II H line is taken as the average of the contaminating flux on both sides of the orders.

All Tables

All Figures

Fig. 1. RV difference between wavelength solutions taken on different days when derived using different algorithms. The RV difference was obtained by comparing all the wavelength solutions with the wavelength solution of reference used in the drift map method, obtained on 2017 November 11. Each column corresponds to a different echelle order, highlighting the behaviour of the wavelength solutions from the blue (order 3) to the red (order 60). Top: RV difference when using only the thorium lines (TH). This is what is done in the old HARPS-N DRS. Middle: RV difference when using thorium lines and the Fabry-Pérot spectrum (TH–FP). Bottom: RV difference when using the drift map method presented in this paper, which means using a reference wavelength solution, and shifting it in velocity after measuring the drift of the thorium lines.

In the text

Fig. 2. RV stability of the wavelength solutions derived with the old (blue) and new (orange) DRS. Top: RV drift of HARPS-N over time, measured by computing the RV drift between all wavelength solutions and the wavelength solution of reference (set for 2017 November 11, see text). We see a general drift over time, interrupted by significant jumps at each warm-up of the detector and power failures (dashed and dotted vertical lines). Bottom: distribution of the RV offsets measured between consecutive wavelength solutions separated by less than 1.5 days, which removes the large offsets observed at each instrument intervention. The MAD of these day-to-day RV offsets is 0.74 and 0.49 m s−1 for the old and new DRS, respectively.

In the text

Fig. 3. Left: number of thorium lines used over time to perform the wavelength solution in the old and new HARPS-N DRS. Here, we show the results for all thorium–thorium and thorium–Fabry-Pérot calibrations. We see that there is a significant variation in the number of lines used in the old HARPS-N DRS, 16%, while in the new HARPS-N DRS, we nearly always used 2041 thorium lines. Right: RV difference between the wavelength solutions obtained with the old HARPS-N DRS and the new one, which uses the drift map method with the new set of thorium lines.

In the text

Fig. 4. Top: variation of the flat fields over time. Here, we show the ratio between the rms and the median error of the residuals obtained by dividing all the flats by one of reference (HARPN.2016-11-20T17:59:35.874_FSPECTRUM_A.fits). Sudden jumps can be explained by periodic detector warm-ups (dashed vertical lines). Bottom: number of raw frames used to build the master flat fields. This number has to be compared to five, which is the number of raw frames included in the old HARPS-N flat field products, and 56, which is the median number of solar observations per day.

In the text

Fig. 5. Raw images from HARPS-N taken before (HARPN.2017-04-09T10-06-46.946, left panel) and after (HARPN.2017-04-12T09-09-40.305, middle panel) a warm-up of the detector, and their ratio (right panel). The location of the core of the Ca II K and H lines in orders one, two, and three is highlighted by grey boxes. As we can see, the core of the Ca II K line in order one (K, bottom left) and the H line in order three (H) are strongly contaminated by ghosts. In addition, we also see that the contamination due to ghosts is significantly reduced with the warm-up process, which removes the humidity accumulated over time on the detector due to a small leak in its cryostat.

In the text

Fig. 6. Top left: comparison between the raw (dots with yellow to blue colour scale) and corrected (red dots) log($R{\prime }_{\mathrm{HK}}$) calcium activity index, as described in the text. Bottom left: differences between the raw and corrected log($R{\prime }_{\mathrm{HK}}$) time series. We clearly see that the slow increase in the activity index followed by sudden jumps at each warm-up of the detector seen in the raw log($R{\prime }_{\mathrm{HK}}$) time series is strongly reduced after correction. Top and bottom right: closer look at the regions around the third and fourth warm-ups.

In the text

Fig. 7. Comparison between the solar RVs reduced using the old HARPS-N (blue) and the new HARPS-N DRS (orange). Top panel: weighted rms for each chunk of RV data of 50 days, which corresponds to two solar rotations and therefore include activity signal. Orange and blue bars correspond to the weighted rms measured on the new RVs and on the old RVs, respectively. The number for each chunk corresponds to the improvement in weighted rms brought by new DRS (in percentage). Second panel: solar RVs as derived by the old and new HARPS-N DRS. The blue and orange lines correspond to a smooth version of the calcium activity index variation linearly fitted to the RVs, to remove the strong long-term correlation between RVs and the log($R{\prime }_{\mathrm{HK}}$), and thus to mitigate the contribution from the solar magnetic cycle (see text). The colour scale for the new RVs corresponds to the S/N of the measurements. We note that we lost a significant amount of data around BJD = 2458100 (December 2017) due to a damaged fibre (see text). Third panel: difference between the new and old solar RVs, for which a long-term trend is clearly visible. We highlighted with red and green vertical lines the significant discrepancies in those RV differences. The red lines correspond to problems in the old DRS RVs, while for the green vertical line, it is not clear if the correct RVs are the ones derived with the old or new DRS. Bottom panel: generalised Lomb-Scargle periodogram of the residuals after removing the long-term trend observed in the RVs. The small inset shows a zoom around the periods 10–30 days to highlight stellar rotation.

In the text

	Fig. 8. Left: histogram of the day-to-day RV offsets measured on the old and new solar RVs after daily-binning the data. Right: intra-day RV variation obtained by removing the daily mean on each single day.
In the text

	Fig. 9. Comparison of the weights for all the lines used to build the CCF between the old and new DRS G2 templates.
In the text

Fig. 10. Top: correlation between the RVs and the log($R{\prime }_{\mathrm{HK}}$) for long-term variations, which in this case means that the time series were smoothed using a 30-day rolling average. Bottom: correlation between the short-term variations, meaning the raw time series minus the smoothed ones. In the legend, we show the values for the R Pearson correlation coefficient for the new and old solar RVs.

In the text

Fig. B.1. Derivation of the wavelength solution for the thorium calibration HARPN.2017-03-25T16-47-52.890.fits using the drift map method with the frame HARPN.2017-11-12T15:39:58.405.fits selected as a reference. Top: first, we measured the drift of all the thorium lines between those two calibrations. Each thorium line is represented as a dot with a size proportional to the error in drift and colour to the drift measured. We note that the thorium lines have been represented at their exact location on the detector using the localisation of the orders (grey dotted lines). Bottom left: then, we fitted a smooth two-dimensional polynomial on the measured drift of all the thorium lines. The corresponding wavelength solution was obtained by shifting the reference wavelength solution with the fitted two-dimensional polynomial. Bottom right: measured drift on the thorium lines between the thorium calibration and the reference before adjusting the two-dimensional polynomial (top) and after fitting for it (bottom). As we see, the χ2 of the two-dimensional polynomial fit is close to one, and the residuals show a scatter, which was measured with the MAD, which is significantly reduced.

In the text

Fig. C.1. Top: thorium flux ratio as a function of time, which is the total flux emitted by the thorium lamp at a given time relative to a reference time. Middle: RV drift of the thorium line at 4833.77 Å measured with the old (blue) and new (orange) DRS. On both data sets, we fitted a model corresponding to a linear drift plus a linear correlation with the thorium flux ratio. In the legend, we report the measured slope of the linear drift. Bottom: RV drift residuals after subtracting the fitted linear drift. In the legend, we report the R Pearson correlation coefficient between those residuals and the thorium flux ratio, as well as the ratio between the rms of those residuals and the median error in RV for this spectral line to detect a non-Gaussian behaviour.

In the text

Fig. C.2. Left: slope measured on the RV drift of all the thorium lines used in the old (blue) and new (orange) DRS. Middle: Pearson correlation coefficient between the RV drift of each thorium line and the thorium flux ratio (see Fig. C.1). Right: rms of the RV drift residuals when subtracting the measured linear drift, divided by the median error in RV. For each subplot, the red vertical dashed lines correspond to the cut-offs used to reject thorium lines showing strong systematics (±10 m s−1 yr−1 for slope, ±0.3 for Pearson correlation coefficient, and < 1.5 for RV residual rms over error in RV). The black histogram corresponds to the selection of lines used 99% of the time, and the orange one corresponds to the final selection used in the new DRS after rejecting lines using the different cut-offs.

In the text

	Fig. C.3. Position of the thorium lines used by the old (blue) and new (orange) DRS as a function of their position on the detector. We note that the position in dispersion direction is in pixels, and the position in the cross-dispersion direction is in echelle order number.
In the text

Fig. D.1. Derivation of the chromospheric emission in the core of the Ca II H line with estimation of the contamination induced by the perturbing ghost. Left: raw image centred on extracted order three with the position of the Ca II H line core. Right: extracted spectra at different positions. The green spectrum is the stellar spectrum extracted from order three (green line on the left panel). The blue and orange spectra are the flux optimally extracted, using the same order profile as order three, minus and plus eight pixels in the cross-dispersion direction on each side of order three (blue and orange lines on the left panel). We note that the extracted fluxes for those spectra have been multiplied by 20 to highlight the contamination from the ghost (red dashed line in the left panel). The vertical green dashed line corresponds to the centre of the Ca II H line core (green dot on the left panel), and the blue and orange vertical dashed lines correspond to the centre of the ghost on both sides of order three (blue and orange dots on the left panel). Assuming that the ghost does not show any curvature, which is a good approximation from one side of the order to the other, we measured the crossing of the ghost and order three (red dot in left panel, and red vertical dashed line in the right panel) as the middle between the orange and blue dots. To measure the chromospheric emission in the core of the Ca II H, the extracted flux from the green spectrum was multiplied with the triangular green response of width 1.09 Å (Vaughan et al. 1978). To measure the effective ghost contamination, we multiplied the contaminating flux on the blue and the orange spectra with the same triangular response, taking into account that the ghost does not cross order three at the centre of the Ca II H line core (the difference is highlighted by the red arrows). Finally, the contamination inside the core of the Ca II H line is taken as the average of the contaminating flux on both sides of the orders.

In the text