A public HARPS radial velocity database corrected for systematic errors

Context. The HARPS spectrograph provides state-of-the-art stellar radial velocity (RV) measurements with a precision down to 1 m/s. The spectra are extracted with a dedicated data-reduction software (DRS) and the RVs are computed by CCF with a numerical mask. Aims. The aim of this study is three-fold: (i) Create easy access to the public HARPS RV data set. (ii) Apply the new public SERVAL pipeline to the spectra, and produce a more precise RV data set. (iii) Check whether the precision of the RVs can be further improved by correcting for small nightly systematic effects. Methods. For each star observed with HARPS, we downloaded the publicly available spectra from the ESO archive and recomputed the RVs with SERVAL. We then computed nightly zero points (NZPs) by averaging the RVs of quiet stars. Results. Analysing the RVs of the most RV-quiet stars, whose RV scatter is<5 m/s, we find that SERVAL RVs are on average more precise than DRS RVs by a few percent. We find three significant systematic effects, whose magnitude is independent of the software used for the RV derivation: (i) stochastic variations with a magnitude of 1 m/s; (ii) long-term variations, with a magnitude of 1 m/s and a typical timescale of a few weeks; and (iii) 20-30 NZPs significantly deviating by a few m/s. In addition, we find small (<1 m/s) but significant intra-night drifts in DRS RVs before the 2015 intervention, and in SERVAL RVs after it. We confirm that the fibre exchange in 2015 caused a discontinuous RV jump, which strongly depends on the spectral type of the observed star: from 14 m/s for late F-type stars, to -3 m/s for M dwarfs. Conclusions. Our NZP-corrected SERVAL RVs can be retrieved from a user-friendly, public database. It provides more than 212 000 RVs for about 3000 stars along with many auxiliary information, NZP corrections, various activity indices, and DRS-CCF products.


Introduction
The High Accuracy Radial velocity Planet Searcher (HARPS, Pepe et al. 2002;Mayor et al. 2003) operates since 2003 at the 3.6 m telescope of the European Southern Observatory (ESO) in La Silla. It is the first fibre-fed high-resolution echelle spectrograph capable of measuring stellar radial velocity (RV) with a precision down to ∼ 1 m s −1 . In this context, HARPS discovered a plethora of exoplanets in the past 15 years. More notably, HARPS proved to be an effective hunter of small and some even potentially temperate exoplanet systems, e.g. around GJ 581 (Bonfils et al. 2005;Udry et al. 2007  . With its unprecedented precision, HARPS is the southern hemisphere backbone Doppler validation instrument for the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), which is uncovering hundreds of small rocky transiting exoplanet candidates around nearby stars. HARPS also offers a large publicly available spectral archive, which has already allowed the post-detection Doppler validation of TESS candidates, for example GJ 143, HD 23472 (Trifonov et al. 2019a), HD 15337 (Gandolfi et al. 2019;Dumusque et al. 2019) and GJ 357 (Luque et al. 2019).
The HARPS spectrograph is precise and stable on years-long timescales thanks to active environmental control (mainly temperature and pressure) and the stability of its atomic standard calibration, typically a ThAr hollow cathode lamp, and since 2011  (Wildi et al. 2010). Although a laser frequency comb is available since April 2015 (Lo Curto et al. 2012) it is not in routine operation yet due to limited spectral coverage and fibre lifetime.
Despite HARPS' stability, its wavelength calibration, drift measurements, and cross-calibration, are nontrivial procedures and are a potential bottleneck limiting the RV precision. The shortcomings in the pipeline, the instrument, or the observations, may lead to systematic errors. An example is the so-called "CCD stitching" systematic, which was discovered with the laser frequency comb (Wilken et al. 2010) and which is not handled by the current pipeline and leads to ∼1 yr periodicity in the HARPS RVs correlated with the barycentric Earth radial velocity (Bauer et al. 2015;Dumusque et al. 2015;Coffinet et al. 2019). Another well-known systematic was introduced with the upgrade of the optical fibres in May 2015 (Lo Curto et al. 2015). This upgrade changed the instrumental profile, and thus the RV offset between the pre-and post-upgrade RVs. This offset is not the same for all stars and might depend on the stellar spectral type.
In this work, we take advantage of the large set of HARPS publicly available wavelength-calibrated spectra accumulated over the years by many programs and groups with different observing strategies and goals. We analyzed the HARPS sample for common RV systematics in an attempt to deliver more precise Doppler measurements for the exoplanet community, and made these data available in a user friendly catalog.
We re-derived more than 212 000 RVs from publicly available stellar HARPS spectra with the SpEctrum Radial Velocity AnaLyser (SERVAL, Zechmeister et al. 2018) pipeline. We then looked for the presence of systematic effects by investigating HARPS' nightly zero point RVs (NZPs): a methodology we had established for CARMENES data ) and for archival HIRES RVs (Tal-Or et al. 2019). Similar approaches had also been applied for the SOPHIE spectrograph (Courcol et al. 2015;Hobson et al. 2018). The systematics revealed in those works can often be traced back and attributed to known instrumental events such as detector changes or fibre coupling problems.
In Sect. 2 we introduce the publicly available HARPS data and the stellar sample we use in our analysis. In Sect. 3 we introduce our HARPS spectral re-processing scheme with the SERVAL pipeline. In Sect. 4 we present our findings regarding HARPS' systematic RV variations. In Sect. 5 we present the main results of our work, and Sect. 6 gives a brief summary and draws some conclusions.

The HARPS data and the stellar sample
Once public, HARPS spectra can be queried and downloaded from the ESO archive using the generic query form 1 , which provides access to all phase 3 data reduced by HARPS' datareduction software (DRS v3.5). The data products also contain detailed auxiliary information regarding the observation such as target's coordinates, Earth's barycentric radial velocity, signalto-noise ratio (S/N) in each order, drift measure, etc. Furthermore, the DRS provides a precise RV measurement derived by the spectrum cross-correlation function (CCF) method using a weighted binary mask (Pepe et al. 2002), as well as the CCF's full-width half-maximum (FWHM) and the Bisector Inverse Slope span (BIS-span) measurements. The FWHM and BIS are often used as stellar activity indicators (e.g. Queloz et al. 2001). 1 http://archive.eso.org/wdb/wdb/adp/phase3_main/form We have developed a new pipeline, which downloads, extracts, and re-processes the large public HARPS archive in a consistent way, and creates an easy user access to its scientific products such as high-precision Doppler measurements and activity indices (see Sect. 5.1). From the ESO-HARPS archive, we have downloaded a total of 264 058 publicly available reduced two-dimensional, multi-order spectra of more than 6 100 objects, which were observed with HARPS between 2003 and mid-2018. We have excluded from our analysis spectra which were not useful for high-precision stellar RV statistical analysis, such as solar spectra, asteroids, the Galilean moons, quasars and supernova candidates, which had been obtained in different astrophysical contexts. Overall, we identified 5260 stellar targets in the HARPS sample. Additionally, we excluded stars with less than three usable spectra -the minimum requirement for creating a meaningful spectral template and thereafter precise RVs with SERVAL. Eventually, we have selected a total of ∼ 3 000 reliable targets of F, G, K, M and L spectral types, to use for our subsequent NZP analysis. Figure 1 shows a Hertzsprung-Russel (HR) diagram of the HARPS sample stars over-plotted on top of the known stars within 100 pc, retrieved from the Gaia DR2 catalog (Gaia Collaboration et al. 2016, 2018. During its operational time, HARPS observed mainly bright nearby main-sequence stars and late-type G and K giant stars. Other stellar objects, such as white dwarfs, very faint late-type M-dwarfs, and brown dwarfs are less suitable for a 3.6 m class telescope and HARPS. Figure 2 describes the HARPS stellar sample in terms of distributions of the estimated stellar distances (Bailer-Jones et al. 2018), V magnitudes and B − V colors. It is evident that the HARPS surveys conducted in the past 15 years focused mainly on nearby main-sequence stars of spectral types G0 V to M6 V, with a median B − V color of 0.722 mag, median apparent V magnitude of 8.4 mag, and with a median distance of ∼120 pc. The stars within 100 pc represent ∼ 67% of our HARPS sample. This collection of stars, are representative of a volume-limited, long-lasting HARPS surveys dedicated to solar-mass G-dwarf stars (Pepe et al. 2004;Lo Curto et al. 2010;Moutou et al. 2011;Lo Curto et al. 2013;Udry et al. 2019), low-mass M-dwarfs Trifon Trifonov et al.: A public HARPS radial velocity database corrected for systematic errors  (Bonfils et al. 2005;Mayor et al. 2009;Forveille et al. 2009;Bonfils et al. 2013;Anglada-Escudé et al. 2016;Astudillo-Defru et al. 2017;Ribas et al. 2018) and metal-poor stars (Santos et al. 2011;Faria et al. 2016;Mortier et al. 2016), which all target nearby stars. The remaining ∼ 28% of the HARPS sample (with distance > 120 pc) are typically bright main sequence stars of spectral types A0 to F6, some fainter and more distant transiting planet hosts observed by more recent RV follow-up campaigns of transit planet candidates from the HATSouth (Bakos et al. 2013;Brahm et al. 2016;Henning et al. 2018;Espinoza et al. 2019), WASP-south (Pollacco et al. 2006;Gillon et al. 2009;Nielsen et al. 2019), and the K2 extended mission (Howell et al. 2014;Grziwa et al. 2016;Johnson et al. 2018), or evolved subgiant and giant branch stars of spectral types G8 IV − K4 III. Figure 3 shows some basic observational statistics from the employed HARPS sample. The left panel shows a histogram of the number of spectra per target, whereas the middle panel shows a histogram of the time baselines of the observed targets. The scatter plot in the rightmost panel shows the relation between the two quantities. It is evident from the Figure that the majority of the targets observed in the sample have fewer than 50 spectra, with about 38% of the targets having fewer than 10 spectra. Some of the targets have very short time baseline below 2 days, but nevertheless often have a sufficient number of spectra. These are likely the result of special observational campaigns such as transit spectroscopy or observations of the Rossiter-McLaughlin (RM;Rossiter 1924;McLaughlin 1924;Queloz et al. 2000) effect of known transiting planets, which require numerous consecutive observations within one or two nights. Most of the stars have time baseline longer than three years and those also tend to have many spectra. The majority of these are most likely known single or multiple-planet hosts and "RV-standard" stars, which show a very low RV scatter, and thus are frequently observed for instrument stability monitoring over the past 15 years of HARPS service. For example, the most frequently observed target with a total of 19 640 spectra that stood out on the top right corner in the right panel of Fig. 3 is our close neighbour α Cen B (Dumusque et al. 2012;Rajpaul et al. 2016). The second most observed star with a total of 11 666 spectra is the RV-standard τ Ceti (HD 10700), which according to Tuomi et al. (2013) has up to five Super-Earth planets. The list of frequently observed stars is followed by another RV standard -Eri (HD 20794) -with a total of 6 775 spectra, which is also a potential multi-planet systems (Pepe et al. 2011;Feng et al. 2017), followed by δ Pav (GJ 780) with 6 003 spectra, and many other long-term observed planetary systems (e.g. Udry et al. 2019). Figure 4 shows statistics of the exposure times and S/N at 550 nm of the public HARPS spectra. The left panel of Fig. 4 shows that in addition to the most commonly used exposure times of 300, 600, 900, 1200 and 1800 seconds, there is an abundance of exposures with less than 150 seconds. Again, these are likely spectra of the most heavily observed RV standard stars, which are very bright, and hence the short exposures and large number of observations. The middle panel of Fig. 4 shows a histogram of the typical S/N at 550 nm calculated by the DRS pipeline. The distribution seems to be bimodal with a peak near S/N ∼ 50, and another broader peak near S/N ∼ 180. Overall, the stars in our HARPS sample are bright (see middle panel in Fig. 2), and the exposure times were selected for achieving the high S/N needed for maximum RV precision. The right panel of Fig. 4 shows a scatter plot of the exposure times and the S/N. The red-dashed lines mark the range of 10 < S/N < 400 adopted in our analysis for the creation of the stellar template needed to derive precise RVs with SERVAL (see Sect. 3).

Deriving RVs and activity indicators with SERVAL
In addition to HARPS DRS RVs, we also derived precise RV measurements with SERVAL ) based on the same DRS spectrum extraction and the same wavelength solution. Instead of using a pre-calculated numeric mask, SER-VAL creates for each observed star a high S/N template spectrum by shifting and co-adding all individual spectra of that star. The template is then used to derive the RVs from the same observed spectra by using a χ 2 -minimization approach. SERVAL is a datadriven approach that aims to exploit all the RV information in a self-consistent way. It post-processes the data, requires at least a few spectra, and provides differential RVs. In contrast, the DRS provides absolute RVs with an excellent precision in an online fashion and given a proper choice of a mask and an initial RV guess. Yet, Anglada-Escudé & Butler (2012) demonstrated that the co-adding method can provide higher RV precision than the CCF method with a weighted binary mask employed in the standard HARPS DRS pipeline, in particular for late-spectral-type stars. We compare the DRS and SERVAL RVs for a subset of the most quiet and heavily observed stars in our sample in Sect. 5.1.
Due to the significant change in HARPS' instrumental profile that accompanied the fibre upgrade in 2015 (Lo Curto et al. 2015), we applied SERVAL separately to spectra obtained before and after the intervention. Hence, for each star that was observed both before and after the intervention, we created two high S/N templates.  . Observation statistics of the HARPS spectra used in this work (∼ 212 000 spectra). The panels show the distribution of the exposure times (le f t), the achieved S/N at 550 nm for the obtained spectra (middle), and a scatter plot of the exposure times versus the S/N (right). The red dashed lines mark the range of S/N adopted in this work to create a template spectrum.
For the template creation and the NZP analysis we use only spectra within the range 10 < S/N < 400. Here we aim to co-add only high-quality spectra and avoid biases from noisy or saturated spectra. We do derive RVs from observations within the range 3 < S/N < 10 and 400 < S/N < 500, but these were flagged as low-S/N and high-S/N observations, respectively, and must be taken with caution. Spectra with S/N < 3 and S/N > 500 were not considered for SERVAL analysis.
The SERVAL pipeline also measures stellar activity indicators, which serve as important diagnostics of the planetary induced Doppler signal hypothesis. For the HARPS spectra SER-VAL provides a measure of emission in the Hα, Na I D, and Na II D activity-related lines. To measure these activity indicators, SERVAL needs an estimate of the stellar absolute RV which is requested from the fits header (here DRS-CCF) or as a fallback from SIMBAD. In addition, SERVAL measures the differential line width (dLW), quantifying variations in the spectral line widths, and the chromatic RV index (CRX) of the spectra, which provides an information about the wavelength dependence of the RV from individual spectral orders as induced by e.g. spots. All the spectral diagnostics obtained with SERVAL come with their uncertainties. For more detailed description of the SERVAL activity indicators we refer to Zechmeister et al. (2018).

Nightly zero-point variations
Since the calibrations for the HARPS observations are typically done at the beginning of each night, we focused our search for systematic effects by calculating the time series of nightly zeropoint RVs (NZPs). To calculate the NZPs we followed a similar procedure to the one described in Tal-Or et al. (2019), where a more detailed description of the algorithm can be found. In short, we calculated a NZP for each night in which at least three different RV-quiet stars (RV scatter < 10 m s −1 ) were observed, by taking the weighted average of the RV, after subtracting from each star its own weighted-average RV (stellar zero point).
Unlike the case of HIRES data, where we had calculated only one NZP time series, for HARPS we had to calculate four different NZP time series: distinguishing between DRS-CCF and SERVAL RVs, and between RVs before the 2015 intervention ("pre RVs") and after the intervention ("post RVs"). The latter separation was also done to prevent the stellar zero-point subtraction from influencing the NZP estimation. Since the intervention introduced spectral-type dependent discontinuous jumps to the stellar RV time series, the stellar zero point of each star in a time series that contains both pre and post RVs heavily depends on the fraction of observations taken before or after the intervention (see also Sect. 4.3). Such inconsistent stellar zero points would add scatter to the NZP time series, which would enhance the NZP uncertainties. In addition, we calculated the four NZP time series for each RV-quiet star individually, after excluding The stellar zero-point subtracted RVs (cyan dots) of all RV-quiet stars have been averaged in each night (NZPs, black points). NZPs calculated with too few RVs (n RV,n < 3, red boxes) or NZPs with too large uncertainties ( 1 m s −1 , red circles) are excluded from correcting the stellar RVs. The RVs in these red-marked nights are corrected by using a smoothed version of the NZP curve, which was calculated with a moving weighted-average (21 d window). For significantly deviating NZPs (magenta circles) we adopted their individual (unsmoothed) NZP value, regardless of their uncertainties. The green dots denote for each night the eventual NZP value that was used for the stellar RV correction. Note the broken x-axis at JD ∼ 2 457 169 (May 26, 2015): since HARPSpre and HARPS-post NZPs were calculated separately, which also included subtracting each star's zero-point RV, both time series are centered around ∼ 0 m s −1 , despite the jump in stellar absolute RVs (see Sect. 4.3). The RV axis was limited to ±10 m s −1 to highlight the small NZP variations. the star itself from the NZP calculation process, to avoid self biasing. Figure 5 shows, separately for SERVAL and DRS RVs, the RVs that were used to calculate the NZPs, the derived individual NZPs, and the NZPs that were actually used to correct the originally-derived RVs. The NZPs that were not used for RV correction are marked in red. These are either NZPs that were calculated with too few RV-quiet stars (n RV,n < 3) or NZPs with too large uncertainties ( 1 m s −1 ). The exact uncertainty threshold that we used to exclude a NZP from the correction stage slightly differs between the four different time series, since it was taken as the scatter of the NZPs around a smoothed version of their time series, which was calculated with a moving weighted-average (21 d window) filter. For correcting the originally-derived RVs, we replaced the excluded NZPs with the filter's value in these nights, and fixed their NZP uncertainties to the scatter mentioned above. On top of that, about two dozen NZPs, significantly deviating from the filter, were adopted for the correction even if their uncertainties were above the threshold. Most of these outlier NZPs come from three week-long observing runs performed in 2013.
This somewhat elaborate NZP correction strategy is tailored to the apparent characteristics of HARPS' systematics. Similarly to Tal-Or et al. (2019), we searched for the typical timescale of HARPS' NZP variations by varying the window size of the filter from 3 to 365 days and looking at the χ 2 dof and the p(F test )value statistics of subtracting the filter from the NZPs. A moving weighted-average filter can be viewed as a non-parametric model of the data, with an effective number of parameters equal to the data time-span divided by the window size. To illustrate the search, Fig. 6 shows the results for HARPS-pre RVs. Unlike the case of the HIRES data, where both statistics indicated a similar typical timescale of 1-2 month (see Fig. 3 Notes. The meaning of the flags: 0 -good NZP (small error and not an outlier); 1 -outlier NZP; 2 -too uncertain NZP, δNZP > max(δNZP); 3 -not enough RV quiet stars to calculate a NZP (n RV,n < 3); 4 -no NZP or filter value were calculated (inside an observing gap > 21 d).
consistent with each other. While the χ 2 dof statistic always preferred a smaller window, pointing towards a stochastic behaviour of the NZPs, the p(F test )-value statistic showed a secondary minimum at a timescale of a few weeks, pointing towards a smooth NZP variation of the instrument. This behaviour probably indicates two effects that drive HARPS' NZP variations: a stochastic one, possibly related to calibration errors, and an additional effect related to more slowly varying instrumental drifts.
The four NZP time series that were actually used to correct the originally-derived RVs are given in an online Table 2 . For each night since JD = 2 452 936 (Oct. 23, 2003), the Table provides its NZP together with the uncertainty (δNZP), the number of RV-quiet stars used to calculate the NZP (n RV,n ), and a flag specifying the type of the NZP. Table 1 summarizes the main characteristics of the four NZP sets. Its upper panel shows the weighted rms (wrms), median NZP uncertainty (med(δNZP)), 2 http://cdsarc.u-strasbg.fr/XXXXXXXXX and the uncertainty threshold that was used to exclude a NZP from being used (max(δNZP)). The lower panel shows the number of nights with a certain flag. The five different flags are explained at the bottom of Table 1. The fact that the NZPs reveal a significant source of systematic RV scatter, on top of the internal uncertainties, is expressed in the ratio wrms(NZP)/med(δNZP), which is in the range of 1.6-2.1 for the four RV sets.

Average intra-night drift
After correcting the RVs for NZP variations, as explained above, we checked the residual RVs for average intra-night drifts. A correlation between the RVs and the time relative to the local midnight t mid is a good indicator not only of actual nightly drifts of the spectrograph, but also of seasonal effects, such as correlations of the RVs with the BERV or the hour angle of observation. Figures 7 and 8 show the four RV-t mid correlations: for preand post-SERVAL RVs, and pre-and post-DRS RVs. The correlations' significance and slopes are given in the insets. While SERVAL RVs show no correlation in pre RVs, and slightly positive correlation in post RVs, the DRS RVs show a small but significant negative correlation in pre RVs, corresponding to an average nightly drift of −1 m s −1 , and no correlation in post RVs.
In many cases the DRS pipeline provides an estimate of the RV drift relative to the nightly calibration sequence, which is calculated by using the light of a reference calibration source injected simultaneously with the star light through a second fibre. As part of deriving the RVs with SERVAL, we corrected the estimated Doppler shifts whenever the drift value was given. The same was done with DRS RVs. Nevertheless, the RV-t mid correlations in pre-DRS and post-SERVAL RVs remain.
We do not know the reason for these correlations, which can be further investigated on a deeper instrumental level. For instance, one can look for correlation with other auxiliary information. For the purpose of this data-driven work, we simply corrected the RVs for the small correlations in all four RV sets. Hence, the final correction model for each RV included both the NZP, according to the night of observation, and the average intranight drift, according to the t mid of observation. The model uncertainties were added in quadrature to the internal RV uncertainties.

The 2015 instrumental RV jump
By the end of May 2015 the HARPS spectrograph underwent a major fibre link upgrade, during which the old circular fibres were replaced by octagonal ones as described in Lo Curto et al. (2015). They find that the main improvement was an increase of the scrambling power by at least a factor of ten, and that due to this any de-centring of the stellar light source on the fibre affect the measured RV by less than 0.5 m s −1 . While this in general is expected to result in more stable RVs, the intervention resulted in a significant change of the instrumental profile. As reported, this introduced an RV offset between the pre-and postupgrade epochs. Lo Curto et al. (2015) investigated this "jump" by comparing the pre-to post-upgrade RVs for more than 20 RV standard stars, and found RV offsets between −2.3 m s −1 and 20.0 m s −1 . Their data also indicates that the offsets might be related to the spectral type of the targets. The reason for the jump is likely the missing optimization of the current DRS extraction pipeline. Maybe it can be better calibrated by a future DRS version, but in particular asymmetry changes in the line spread func-Article number, page 6 of 17   tion are challenging to model. A dependence on spectral type is expected, since the systematics might vary across the detector, and the RV information also changes with spectral type across the detector. We also investigated the magnitude of the jump with SER-VAL. For this purpose, we recomputed again the HARPS RVs, but this time with a common stellar template, meaning that preand post-data were processed jointly and not separated. Then we recomputed the NZPs for this entire data set as described in Sect. 4.1.
Following this procedure, we found that while the computed NZPs were stable for the pre-upgrade epoch, they showed a steady drift for the post-upgrade period. Therefore, in order to estimate the RV offset, we applied separate linear fits to the NZPs before and after the intervention break. Only valid NZPs that were computed from nights with sufficient RV data were used for the fits and consequently for the determination of the offset. This procedure, using the entire data set, yielded a jump of 8.19 ± 0.38 m s −1 between the last pre-upgrade night (JD = 2 457 163) and the first post-upgrade night (JD = 2 457 173) 3 .
In order to minimize systematics in the offset estimation that might be introduced by any drifts of the NZPs, in particular in combination with unequal distributions in the numbers of observations between the two epochs, we repeated the same procedure (NZP computation and jump determination) with a modified (reduced) data set. This "stripped" set was created by discarding for each star some of the earliest or latest RVs, so that the numbers of nights in which the stars were observed, were equal for the two epochs. The NZPs estimated from this stripped data are illustrated in Fig. 9. The plot is clipped to a time span of four years around the instrumental intervention. The jump in the NZPs between nights prior to the upgrade and those after is evident and is also reflected in the bimodality of the NZP distribution. For the stripped data we found an RV offset of 7.93 ± 0.33 m s −1 , which is consistent with the first approach. The reported errors are propagated from the linear fits, where the uncertainties of the regressions were estimated by bootstrapping, for which the data was re-sampled 1000 times. For comparison we repeated the entire analysis also with the DRS RVs. With offsets of 10.21 ± 0.43 m s −1 for the entire data set, and 9.35 ± 0.38 m s −1 for the stripped data set, they appear to show higher jumps than those determined from the SERVAL RVs.
We also investigated the dependence of the magnitude of the RV jump on spectral type. To do so, we gropued the stars by Article number, page 7 of 17 A&A proofs: manuscript no. Trifonov_NZP_HARPS  their spectral type, and repeated the analysis on each sub-sample alone. The results from all analyses are summarized in Table 2.
For the sub-samples of later spectral types the offset estimates from both methods (full data and stripped data) are similar and consistent within their uncertainties, while they deviate a bit for spectral type G. No jump values could be estimated for the stripped data sets of the F-type stars. This is due to the combination of a small number of stars and the reduced data set, which together lead to insufficient number of nights for which the NZPs could be estimated. However, they were determined for the full data sets. The dependence on the spectral type is evident. The offsets are highest for early type stars and appear to decrease for lower effective temperatures. For M dwarfs the sign of the jump is inverted. The apparent relationship is illustrated in The resulting offsets should be treated with care, as it is expected that the jumps vary from one source to another. However, the results presented here are still a valuable source of information. They can either be used for sanity checks, when offsets between pre-and post-upgrade RV data are treated as free parameters in RV modeling, as usually done in planet searches, or they can even serve as priors for the offset. In particular for large amplitude signals, such priors can serve as useful constraints.

The new HARPS RV database
Butler et al. (2017) published an example of a well documented, user-friendly, high-precision RV database. Their database contains a large collection of precise Doppler measurements, derived from ∼ 65 000 spectra of ∼ 1700 stars, as well as stellarline activity-index measurements, obtained with the iodine cell method (Marcy & Butler 1992;Valenti et al. 1995;Butler et al. 1996) between 1996 and 2014 with the KECK-HIRES 4 spectrograph (Vogt et al. 1994). The RVs published by Butler et al. (2017) were by far the most precise and extensive HIRES-RV archive available to the exoplanet community. These public HIRES RVs were the basis for a number of new exoplanet discoveries and orbital updates (e.g. Butler et al. 2017;Trifonov et al. 2018;Kaminski et al. 2018;Trifonov et al. 2019b;Tuomi et al. 2019), and form an important RV validation archive for the TESS candidates in the Northern hemisphere. The large sample size of the HIRES RV archive has also allowed us to identify and correct the data for systematic variations, by calculating nightly zero point variations, and an average intra-night drift (Tal-Or et al. 2019). Following a similar route to Butler et al. (2017), we created HARPS-RVBank -a public database based on the results presented in this work. HARPS-RVBank is available on CDS 2 or on its official web page 5 and provides up-to date SERVAL and DRS data products for the HARPS targets. Table A.1 shows an excerpt of the database with some important columns. The HARPS-RVBank provides original as well NZP-corrected SER-VAL and DRS RV measurements, their BJD epoch, activity index measurements, and uncertainties. From the DRS products, the user can find the CCFs FWHM, contrast, and the BIS-span measurements. From the SERVAL spectral analysis, we provide the chromatic index (CRX), the differential line width (dLW), and the Hα, Na I D, and Na II D activity-related line measure-ments. We note that the activity time series are also affected by the 2015 fibre upgrade (see Sect. 4.3) leading to a notable postupgrade "jump" in the data. We recommend to treat the pre-and the post-upgrade activity time series as taken from two different instruments, before testing the data for the presence of periodic signals.
The new archive also provides all applied individual RV corrections to the data, such as barycentric Earth radial velocity (BERV), secular acceleration of the RV (SA; Kürster et al. 2003), Fabry-Perot (FP) drift, DRS and SERVAL NZP time series, and the final correction value of each RV, including the average intranight drift. Additionally, we provide for each epoch auxiliary observational information such as exposure time, S/N of the spectra at 550 nm, quality flag, type of DRS binary mask used, principal investigator (PI) and ESO program-ID of the observation.
For user flexibility we also provide a github repository 6 of the HARPS-RVBank database, where the user can find the final products, useful tools, and instructions. There we also provide a MySQL, and a JavaScript Object Notation (JSON) versions of the HARPS-RVBank database, which could be easily integrated as an object storage in modern programming languages such as Python, or used in web interfaces. Figure 11 compares SERVAL and DRS-CCF as RV derivation tools. For the comparison, we focus on the most quiet stars (RV scatter < 5 m s −1 ), and compute the average (and median) ratio std(RV-SERVAL)/std(RV-DRS), where we use wrms as our std estimator. We find that for HARPS-pre SERVAL yields RVs with an average wrms improvement of ∼ 5% (∼ 4% median), while for HARPS-post the average wrms improvement is ∼ 4% (∼ 2% median). Hence, for most stars SERVAL yields a slightly better RV precision, and should in general be preferred over the nominal DRS RVs. Figure 12 shows the impact of the NZP correction on the wrms of all RV-quiet stars (RV scatter < 10 m s −1 ). Comparing the impact of the correction between HARPS-pre and HARPSpost, we see that the average impact on HARPS-pre data is smaller than the average impact on HARPS-post data, for both SERVAL and DRS RVs. In order to quantify the effect, we again look at the average (and median) wrms improvement due to the NZP correction. For both SERVAL and DRS RVs, the NZP correction yields only a marginal effect on HARPS-pre data (< 2%), while for HARPS-post data the average reduction of wrms is ∼ 8% (∼ 6% median). This finding is consistent with the fact that, for both RV-derivation pipelines, HARPS-post NZPs have larger scatter and smaller uncertainties than HARPS-pre NZPs (see Table 1), leading to a more significant correction of post RVs with the calculated NZPs.

DRS versus SERVAL and the impact of the NZP correction on RV-quiet stars
Interestingly, in HARPS-pre there is a larger number of stars with a correction that significantly deviates from the mean. Specifically, in HARPS-pre there are ∼ 20 stars with a wrms improvement of 2 m s −1 , while there are only four such stars in HARPS-post. This effect is probably related to the larger number of significantly outlying NZPs in HARPS-pre (see Table 1). For instance, there are two stars in HARPS-pre with a wrms reduction from > 5 m s −1 to < 2 m s −1 : HD 145927 and HD 197818. HD 145927 has 10/37 of its RVs corrected by outlier NZPs of ∼ −5 m s −1 , and HD 197818 has one RV corrected  Table 1.  Figure 13 shows the impact of the NZP correction on the median RV uncertainty per star (med(δRV)). Co-adding the NZP uncertainties (δNZPs) to the original RV uncertainties enlarged the med(δRV) values. However, the estimated δNZPs rarely exceeded ∼ 1 m s −1 . We believe that the new RV uncertainties better represent HARPS' RV precision, and will require a smaller jitter term in modeling the RV time series with Keplerian orbits. Moreover, the nights with the largest number of bright RV-quiet stars, observed under good conditions, naturally give the smallest δNZPs. This, in turn, gives higher weight to the RVs from the nights in which we have a good estimate of the calibration errors.
Another demonstration of the importance of the NZP correction can be seen in Fig. 14, which compares the RV wrms per star before and after the 2015 intervention. Without NZP correction, HARPS' fibre exchange actually worsened the wrms per star by ∼ 10% on average. However, after correcting for the NZPs, we find it to improve the wrms per star by ∼ 10% on average. Moreover, the pre and post wrms values are more correlated when correcting for the NZPs. This indicates that the fibre exchange indeed improved the instrument performance, but the wrms values of the most quiet stars were still dominated by instrumental systematic effects. It is the NZP correction that enables the better precision achieved with HARPS' new fibres.
Focusing again on the most quiet stars (RV scatter < 5 m s −1 ), Fig. 15 shows end-to-end comparison between the RV performance of the nominal DRS-CCF (without NZP correction) and SERVAL with NZP correction. For HARPS-pre data, the average wrms improvement is ∼ 5% (∼ 4% median), which is dominated by a small number of stars with a relatively large wrms improvement. For HARPS-post data, the average wrms improvement is ∼ 15% (∼ 18% median), with only a few stars above the 1:1 line. The difference can be explained by the different NZP behaviors of HARPS-pre and HARPS-post RVs. We conclude that for HARPS spectra the NZP-corrected SERVAL RVs are in general more precise than the DRS RVs, and we regard them as the main product of this work.

A practical example: Orbital update of the GJ 253 multi-planet system
We now make a direct comparison of the official HARPS-DRS RV data, and our final product -the SERVAL NZP-corrected RV data, by testing their overall quality when modeling an actual planetary system discovered with HARPS. For this purpose we use the known multi-planet system GJ 253 (HD 51608), whose RVs are consistent with two Neptune-mass planets Udry et al. 2019). We have selected this system rather arbitrarily, as in our HARPS archive we find many known planetary systems, which could serve as better (or worse) examples for testing our SERVAL-NZP data. The GJ 253 system, however, is appropriate in the context of the results presented in this work, because: -The DRS data have no outliers (e.g. wrong user RV-guess), meaning that we can perform a full data set comparison with the SERVAL data. -Data were taken both before (218 RVs) and after (9 RVs) the May 2015 fibre upgrade. -GJ 253 has a sufficient number of archival RV data for a proper statistical comparison. -The planetary signals are strong and unambiguous (i.e. likely real, and not due to stellar activity). -The system is fairly complex (a two-planet system), but is still easy to analyze, since there is no need of N-body modeling. -We can perform an update of the systems' orbital solution. There is no evidence of significant GLS power in the activity indices from DRS and SERVAL, that could be associated with the two RV signals. The only exception is the SERVAL's NaD I activity indicator, which seems to show a marginally significant power near 14.04 d. This NaD I peak, however, is only the tenth strongest peak in the NaD I GLS power spectrum. A closer inspection of the RV and the NaD I time series, showed that these two peaks are sufficiently distant (>3σ) in frequency space, and show no correlation. Therefore the strong RV signal, and the marginally significant NaD I signal are likely not related to each other.
For the fitting analysis of the GJ 253 system we use The Exo-Striker fitting toolbox 7 (Trifonov 2019). To identify the planetary signals embedded in the data, we perform a "blind-search" using the The Exo-Striker's "Auto fit" algorithm, which, as its name suggests, automatically inspects the RV data for periodic signals via GLS periodogram search and performs pre-whitening signal subtraction (Hatzes 2013). Finally, when no significant peaks are left in the RV data residual periodogram The Exo-Striker performs a subsequent simultaneous best-fit optimization of the planetary semi-amplitude K, the orbital period P, eccentricity e, argument of periastron ω, and mean anomaly M at the first observational epoch. We consider the HARPS-pre and HARPS-post fibre upgrade data as obtained from separate instruments, and thus we optimize their RV offsets and RV velocity jitter (Baluev 2009) simultaneously with the planetary parameters. Using this approach we were able to instantly identify the planetary signals published in the discovery paper by Mayor et al. (2011), in both the DRS data and the SERVAL NZP-corrected data. As a next step, The Exo-Striker employs a Markov Chain Monte Carlo (MCMC, via the emcee sampler; Foreman-Mackey et al. 2013) sampling around the best fit, and we adopt the 1σ confidence level of the posterior distributions as parameters uncertainties. Table 3 summarizes the best-fit parameters and MCMC derived asymmetric uncertainties from our dual RV analysis, while Fig. 17 shows the best two-planet fits of the DRS and SERVAL NZP-corrected RV data sets. Both data sets lead to consistent best-fit planetary periods of P b ∼ 14.07 days and P c ∼ 95.9 d, which are in agreement with the period estimates by Mayor et al. (2011), who gives P b = 14.070±0.004 days and P c = 95.42±0.39 days. We derive lower and better constrained eccentricities when compared to those provided in the discovery work of Mayor et al. (2011). The DRS data leads to e b = 0.087 +0.030 −0.056 and e c = 0.258 +0.059 −0.099 , SERVAL NZP-corrected data suggest e b = 0.091 +0.027 −0.061 and e c = 0.227 +0.050 −0.074 , while Mayor et al. (2011) gives e b = 0.15±0.06 and e c = 0.41±0.18. We note that Mayor et al. (2011) had used only a little more than half of the HARPS data we use in our analysis, and thus it is not surprising that our estimates are better constrained by the larger, and longer baseline data set. Udry et al. (2019) have recently published an orbital update of the GJ 253 system. They based their orbital solution and activity analyses only on HARPS spectra obtained before the fibre upgrade, but with data reduction done using a new version of the HARPS-DRS pipeline. The MCMC-based best-fit solution given in Udry et al. (2019) is generally consistent with our bestfit estimates (see Table 10 in their paper).
To derive the planetary minimum masses and semi-major axes, we adopt the recent stellar mass estimates for GJ 253 by Soto & Jenkins (2018), who used HARPS spectra to infer stellar parameters and derived a stellar mass of M = 0.85 M . The DRS data and the SERVAL-NZP corrected data mutually agree on the RV semi-amplitudes of GJ 253 b and c, and thereafter on the derived minimum masses of m b sin i ∼ 0.04M jup and m c sin i ∼ 0.05M jup . Thus, our orbital update analysis is consistent with two warm Neptune-mass planets in orbit around GJ 253, as it was reported by Mayor et al. (2009), andUdry et al. (2019).
The orbital dynamics of the GJ 253 system is beyond the scope of this work, but for the sake of completeness we also examine the long-term orbital stability of the best-fit configurations achieved form the different data sets. For this purpose, we adopt the Wisdom-Holman (also known as MVS; Wisdom & Holman 1991) N-body integrator implemented in the The Exo-Striker, which includes General Relativistic (GR) precession correction term. We find that given the small orbital separations and significant planetary eccentricities observed in the GJ 253 system, the GR precession effects on the orbital dynamics are significant, and thus the GR correction must be included in the long-term evolution of the system to assure realistic dynamical outcome. Overall, a crude stability analysis shows that the GJ 253 system is stable for at least 10 M yr, exhibiting an interesting apsidal alignment libration ∆ = bc ∼ 0 • with a semi-amplitude of 60 • , and significant oscillations of the planetary eccentricities in the range 0.02 < e b < 0.17 and 0.21 < e c < 0.23 with mean values of e b = 0.11, and e c = 0.22, respectively. The time scale of the secular orbital oscilation is ≈ 37 000 yr. In Fig. 16 it is evident that the SERVAL dLW, CRX, and H α activity indicators suggest a periodicity near 36.5 d. It is particularly strong in the SERVAL dLW, where this periodicity has a significant power. A peak near this period is also detected in the DRS BIS-span, and the DRS FWHM time series. A peak at a similar frequency was also detected by Udry et al. (2019) on their HARPS activity time series. We are likely witnessing stellar spots rotating with the rotational period of the star. Interestingly, the RV residuals in both DRS and SERVAL show a strong, but insignificant GLS peak near 38.9 d. While this peak seems sufficiently distant from the ∼ 36.5-d period seen in the activity index Article number, page 12 of 17 Trifon Trifonov et al.: A public HARPS radial velocity database corrected for systematic errors time series, it is still possible that it could be related to the stellar rotation. Finally, when it comes to the quality of the two fits, the DRS-CCF data yields ∆ ln L = 180.00, with respect to the null hypothesis (i.e. no planets), while the SERVAL-NZP has ∆ ln L = 185.99, which means that the SERVAL-NZP corrected RV data adds significant evidence in favour of the two planets. This also manifests in a lower weighted rms of the SERVAL-NZP residuals compared to DRS-CCF (see Table 3). Therefore, our orbital update on the GJ 253 system is based on the two-planet Keplerian modeling of the SERVAL NZP-corrected data.
This practical example of the GJ 253 system shows that our SERVAL data are indeed a better choice with respect to the official DRS data. It is still possible that for other systems the DRS data would lead to better fits, but overall, given the sample statistics comparison given in Sect. 5.2 (see also Figs. 11 and 15), we are confident that the NZP-corrected SERVAL data should be preferred in most cases. In particular for stars of late spectral type such as K, M and L, we expect SERVAL to outperform the DRS-CCF.

Summary and conclusions
In this paper we present an independent systematic analysis of the HARPS spectral archive. In particular, we re-calculated Doppler velocity measurements of the publicly accessible spectra, and performed stellar line analysis, which is important for validating planetary induced Doppler signals. For this part of the analysis, we applied the SERVAL RV pipeline, with which we derived slightly more precise RV measurements when compared to those derived by the official ESO-DRS pipeline. We find that for a sub-sample of stars with a very small RV scatter (< 5 m s −1 ), SERVAL RVs are more precise than DRS RVs by ∼ 5%, on average.
We make all of our HARPS results publicly available, as a service to the exoplanet community. We provide original uncorrected DRS and SERVAL pipeline data, and self-corrected (for NZPs and average intra-night drifts) DRS and SERVAL RVs.
All relevant results of our study are made public in the user friendly database "HARPS-RVBank", which is available on the cds, github, or as a stand alone webpage where the user can browse for data.
This makes our HARPS-RVBank the first easy-to-access publicly available HARPS RV archive, which to our knowledge contains the most precise RV data products to date.
Another objective of our work was to study whether HARPS Doppler measurements suffer from systematic errors, which could bias the orbital parameter estimates of small planets, or worse, induce spurious planetary discoveries. We find that despite being, with no doubt, a state-of-the-art RV instrument, HARPS also suffers from small but significant systematic errors of the instrumental zero-point. Our NZP analysis reveals stochastic zero-point variations of ∼ 1 m s −1 , smooth zero-point variations, with a magnitude of ∼ 1 m s −1 and a typical timescale of few weeks, and a few dozen nights whose NZPs significantly deviate from the general zero-point trend by few m s −1 (and up to 30 m s −1 ). In addition, we find small ( 1 m s −1 ) but significant intra-night drifts in DRS RVs before the 2015 intervention and in SERVAL RVs after it. The HARPS NZPs systematic errors are likely related to the non-stability of the daily wavelength calibration (Dumusque 2018). The next DRS version will likely have an improved wavelength calibration (Coffinet et al. 2019), which to a large degree might resolve these systematic errors. Then, it would be interesting to recompute the NZPs as a quality check of the new DRS wavelength calibration scheme.
Correcting HARPS RVs for the systematic's model, we find additional improvement of the RV scatter, mainly for observations after the 2015 intervention. Considering the combined effect of deriving the RVs with SERVAL, and correcting them for the small systematic errors, we find the average wrms improvement of pre RVs to be ∼ 5%, and of post RVs to be ∼ 15%. For a small number of stars, whose observations were most affected by the significantly-deviating NZPs, we find a much more significant wrms improvement, by a factor 2.
Investigating the RVs of a sub-sample of RV-quiet stars that were observed both before and after the 2015 HARPS optical fi-A&A proofs: manuscript no. Trifonov_NZP_HARPS bre upgrade, we find a discontinuous jump in their absolute RVs that is independent of the RV derivation software (i.e. we find similar results with DRS and SERVAL). Similarly to Lo Curto et al. (2015), we find the jump to be strongly dependent on the spectral type of the target: from ∼ 14 m s −1 for late F-type stars, to ∼ −3 m s −1 for late M dwarfs. As a demonstration of the new data quality, we provide new orbital estimates of the GJ 253 multi-planet system based on our new HARPS-SERVAL NZP-corrected data, updating the planetary minimum masses and orbital elements. Similarily to Udry et al. (2019), we show that the GJ 253 b & c orbits are probably less eccentric than was previously estimated, when fewer RV data were available. This shows that it is important to update the orbital elements of known planetary systems when more data are accumulated, as this might remove possible higher-eccentricity biases, for example. This is especially valid for multi-planet systems. In those systems the eccentricity is an important parameter that can determine their dynamical properties and might shed some light on their formation and evolution.
The HARPS-RVBank is a valuable data source for planet search, re-analysis of known planetary systems, and validation of newly discovered transiting planets. For its better precision, we strongly recommend using the NZP corrected SERVAL RVs, which is the main product of this work. seminal work and consistent effort, this study could not have been performed. We are also extremely grateful to all the PIs and observers associated with the ESO programmes listed below. We respect their hard work, and we hope that they will benefit from our recomputed data. T.T. thanks to Thomas Henning and Martin Kürster for support during this study. L.T. and S.Z. acknowledge support from the Israel Science Foundation (grant no. 848/16). M.Z. is supported by the Deutsche Forschungsgemeinschaft under DFG RE 1664/12-1 and Research Unit FOR2544 "Blue Planets around Red Stars", project no. RE 1664/14-1. T.T. acknowledges support by Bulgarian National Science Programme "Young Scientists and Postdoctoral Candidates 2019". This work is based on observations collected at the European Organization for Astronomical Research in the Southern Hemisphere under ESO programmes: 0100.C-0097, 0100.C-0111, 0100.C-0414, 0100.C-0474, 0100.C-0487, 0100.C-0750, 0100.C-0808, 0100.C-0836, 0100.C-0847, 0100.C-0884, 0100.C-0888, 0100.D-0444, 0100.D-0717, 0101.C-0232, RV (e) σRV (