© ESO, 2017

1. Introduction

The radial velocity (RV) method has been a very successful technique to discover and characterise stellar companions and exoplanets. There are many algorithms to compute RVs, which can be grouped depending on their complexity and choices for the model of the reference spectra, number of model parameters, and statistics; these algorithms range from simple cross-correlation with binary masks (Queloz 1995; Pepe et al. 2002), to least-squares fit with coadded templates (Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015), and finally to least-squares fitting with modelling of line spread functions (Butler et al. 1996; Endl et al. 2000). Recently, Gaussian processes have also been proposed to derive RVs (Czekala et al. 2017). The algorithm choice is influenced by instrument type (e.g. those with stable point spread functions) and the need for accuracy (e.g. fast and flexible reduction directly at the telescope, generality) or high precision (David et al. 2014). While most algorithms compute radial velocities in the wavelength or pixel domain, there are also methods that work in the Fourier domain (e.g. Chelli 2000) and can efficiently disentangle double-lined spectroscopic binaries (KOREL; Hadrava 2004). A review can be found in Hensberge & Pavlovski (2007).

In this work we present our methods to measure RVs and spectral diagnostics, which are implemented in a Python programme called SERVAL (SpEctrum Radial Velocity AnaLyser). It aims for highest precision with stabilised spectrographs. Therefore, we employ forward modelling in pixel space to properly weight pixel errors, and we reconstruct the stellar templates from the observations themselves to make optimal use of the RV information inherent in the stellar spectra (Zucker 2003; Anglada-Escudé & Butler 2012).

The SERVAL code was developed as the standard RV pipeline for CARMENES, which consists of two high resolution spectrographs covering the visible and near-infrared wavelength ranges from 0.52 to 1.71 μm (Quirrenbach et al. 2016). The CARMENES consortium regularly obtains spectra for about 300 M dwarfs to search for exoplanets. The spectra are processed, wavelength calibrated, and extracted by a pipeline called CARACAL (Caballero et al. 2016). Using CARMENES and HARPS data we validate the performance of SERVAL, which is publicly available ¹.

2. Radial velocity algorithm

Our RV computation aims for highest precision and is based on least-squares fitting. Anglada-Escudé & Butler (2012) demonstrated that this approach can yield more precise RVs than the cross-correlation function (CCF) method. This is because the RV precision depends on both the signal-to-noise ratio of the data and the match of the model to the data. With a proper model of the observations, the RV information, i.e. each local gradient, is optimally weighted (Bouchy et al. 2001) and outliers can be detected.

Our very general concept is to decompose simultaneously all observations into (1) a high signal-to-noise template F(λ); (2) RV shifts v_n for observation number n; and (3) multiplicative background polynomials p(λ) to account for flux variations, e.g. wavelength-dependent fiber coupling losses due to imperfect correction of the atmospheric dispersion. Therefore, our forward model f(λ) for the observed flux is $\mathit{f}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathit{p}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{·}\mathit{F}\mathrm{\left(}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\right(}\mathit{\lambda ,v}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(1)where the Doppler equation λ′(λ,v) is given in Eq. (6).

We denote some n = 1... N observations taken at times t_n each consisting of flux measurements f_n,i at pixel i with calibrated wavelengths λ_n,i and flux error estimates ϵ_n,i. For cross-dispersed echelle spectrographs the measurements are carried out in several echelle orders o (so a more explicit indexing is f_n,o,i). The weighted sum of residuals is $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}& \mathrm{=}& \\ & \mathrm{=}& \end{array}$with weights ${\mathit{w}}_{\mathit{n,i}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{ϵ}}_{\mathit{n,i}}^{\mathrm{2}}}$, the polynomial coefficients a, and the coefficients b describing the template.

However, this simultaneous approach is hardly feasible in practice because of the large amount of data (spectra, orders, pixels, N × O × K) and large number of parameters. Therefore, we perform the decomposition sequentially and iteratively. First, the observed spectrum with the highest signal-to-noise is taken as a reference to shift all observations into this reference frame and to coadd them to a high signal-to-noise template F (Sect. 2.3). With this new template the radial velocities are recomputed while the coefficients of the background polynomial are fitted simultaneously (Sect. 2.2). This one iteration is usually sufficient (Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015). Before we detail the two steps in Sects. 2.2 and 2.3, we define the model and prepare the data.

2.1. Model details, data, and input preparation

Given a template F(λ′) of the emitting source, a spectral feature at reference wavelength λ′ in F appears at λ in the observation f, i.e. f(λ) = F(λ′). The measured wavelength λ depends on the radial velocity of the moving source, i.e. the star, and on the velocity of the moving observer on Earth (Wright & Eastman 2014).

To eliminate the contribution from Earth’s motion, we transform the measured wavelengths on Earth to the solar system barycentre via the barycentric correction (cf. Eq. (25) in with λ_meas = λ and λ_meas,true = λ_B) $\begin{array}{ccc}{\mathit{\lambda }}_{\mathit{B}}\mathrm{=}\mathit{\lambda }\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathit{B}}\mathrm{)}\mathit{,}& & \end{array}$(4)where z_B is the total redshift due to barycentric motion of the observer. To compute the barycentric correction, SERVAL uses the Fortran based code of Hrudková (2009)², which we modified to account for the proper motion of the star. Moreover, we correct for secular acceleration (Zechmeister et al. 2009), which is in particular important for high proper motion stars and requires the parallax of the star. For ultimate precision in the barycentric correction, the code from Wright & Eastman (2014) with 1 cm/s accuracy can be used instead, which includes relativistic corrections ³.

The stellar radial motion v = cz redshifts the wavelength (cf. Eq. (7) in Wright & Eastman 2014) ${\mathit{\lambda }}_{\mathit{B}}\mathrm{=}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{,}$(5)where c is the speed of light. Therefore, the equation to Doppler shift the template with radial velocity parameter v is $\begin{array}{ccc}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\lambda ,v}\mathrm{\right)}\mathrm{=}\mathit{\lambda }\frac{\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathit{B}}}{\mathrm{1}\mathrm{+}\frac{\mathit{v}}{\mathit{c}}}\mathrm{·}& & \end{array}$(6)With the approximation $\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{\mathit{v}}{\mathit{c}}}\mathrm{\approx }\mathrm{1}\mathrm{-}\frac{\mathit{v}}{\mathit{c}}$ and barycentric corrected wavelengths, one could write this as ${\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{=}{\mathit{\lambda }}_{\mathit{B}}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}\frac{\mathit{v}}{\mathit{c}}\mathrm{)}$ as, for example, in Anglada-Escudé & Butler (2012). The approximation is accurate to 1 m/s over a velocity range of 17.3 km s^-1, i.e. well sufficient for the small amplitudes of exoplanets. But since there is actually no need for this approximation, we employ the exact equation.

The SERVAL code has the option to perform the computation in logarithmic wavelengths. Then Eq. (5) becomes $\ln {\mathit{\lambda }}_{\mathit{B}}\mathrm{=}\ln {\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{+}\ln \left(\mathrm{1}\mathrm{+}\frac{\mathit{v}}{\mathit{c}}\right)\mathrm{·}$(7)For small shifts the Doppler shift is linear, $\ln {\mathit{\lambda }}_{\mathit{B}}\mathrm{\approx }\ln {\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{+}\frac{\mathit{v}}{\mathit{c}}$. This approximation is accurate to 1 m/s over a velocity range of 24.5 km s^-1 and should be avoided in high RV precision work in particular in the barycentric correction in Eq. (4).

We assume that the spectrum is continuous and use cubic spline interpolation to evaluate the template F at any wavelength λ, which is needed for the forward modelling in Eq. (1).

Furthermore, we generate and propagate a bad pixel map to flag pixels with saturation, significant negative flux (f_i,n< −3ϵ_i,n), significant deviation in the fitting (outliers), or contamination by tellurics and sky emission lines. The default telluric mask flags atmospheric features deeper than 5%. We also check the fits header for the correct observing mode (e.g. accuracy and efficiency mode in HARPS, dark time), available drift measurements, and suitable signal-to-noise ratios. Those cases are flagged and, when required, excluded in the analysis.

2.2. Least-squares RVs

For the determination of the RV, we optimise Eq. (3) with respect to the polynomial coefficients a and the RV shift v. The model is linear in a, but the RV parameter v makes the fit non-linear. There are several algorithms, e.g. Levenberg-Marquardt or bisection, to solve non-linear least-squares problems. We use, in analogy to the CCF method, a direct stepping through the velocity parameter space v with a default step size of Δv = 100  m/s as a compromise between oversampling the resolution element (~km s^-1) and computational speed. At fixed velocity v_k, the Doppler-shifted template F_i,k = F(λ_i,v_k,b) is evaluated at each pixel i, and χ²(v_k) is obtained from a simple linear least-squares fit for a${\mathit{\chi }}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}}{\mathit{w}}_{\mathit{i}}\mathrm{[}{\mathit{f}}_{\mathit{i}}\mathrm{-}\mathit{p}\mathrm{(}{\mathit{\lambda }}_{\mathit{i}}\mathit{,a}\mathrm{)}\mathrm{·}{\mathit{F}}_{\mathit{i,k}}{\mathrm{]}}^{\mathrm{2}}\mathit{.}$(8)In the form ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}{\sum }_{\mathit{i}}{\mathit{F}}_{\mathit{i}}^{\mathrm{2}}{\mathit{w}}_{\mathit{i}}\mathrm{[}\frac{{\mathit{f}}_{\mathit{i}}}{{\mathit{F}}_{\mathit{i}}}\mathrm{-}\mathit{p}\mathrm{(}{\mathit{\lambda }}_{\mathit{i}}\mathit{,a}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}$ and with the substitution ${\mathit{F}}_{\mathit{i}}^{\mathrm{2}}{\mathit{w}}_{\mathit{i}}\mathrm{\to }{\mathit{w}}_{\mathit{i}}$ and $\frac{{\mathit{f}}_{\mathit{i}}}{{\mathit{F}}_{\mathit{i}}}\mathrm{\to }{\mathit{f}}_{\mathit{i}}$, standard library routines for polynomial fitting can be applied (division by zero flux F_i needs to be handled). We include the polynomial in the model. Anglada-Escudé & Butler (2012) argued that this “would couple the flux normalisation coefficients to Doppler factor in a non-linear fashion” in their least-squares algorithm with iterative differential corrections and, therefore, they applied it to the data. However, they do not correspondingly re-adjust the errors as implied in our derivation here. Since in practice the variations in the polynomial function should be only a few percent, their approximation is feasible.

Fig. 1 Least-squares RV fit for RV with CARMENES VIS data for GJ 1265. Top: one observation (red) and the best fit model (blue) are shown. Data points within telluric regions (grey area) are a priori excluded. Data points classified as outliers are indicated in cyan. Bottom: fit residuals are shown.

Fig. 2 Spectra coadding of seven CARMENES VIS spectra of GJ 1265. The seven observations (colour coded), each normalised by a background polynomial, are fitted by a uniform cubic B-spline (black) where basically each knot value is a free parameter (black dots). Points affected by tellurics or outlying by 5σ are excluded from the fit (open circles).

The χ²(v) function sampled at $\mathrm{(}{\mathit{v}}_{\mathit{k}}\mathit{,}{\mathit{\chi }}_{\mathit{k}}^{\mathrm{2}}\mathrm{)}$ is explored for its global minimum. Around its minimum, the χ² function takes a parabolic shape and the first derivative vanishes, (χ²)′ = 0. The Taylor expansion of the χ² function at the minimum to second order is ${\mathit{\chi }}^{\mathrm{2}}\mathrm{(}\mathit{v}\mathrm{+}\mathit{\delta v}\mathrm{)}\mathrm{\approx }{\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{\chi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{\prime }}\mathit{\delta v}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{\chi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{\prime \prime }}\mathit{\delta }{\mathit{v}}^{\mathrm{2}}\mathrm{=}{\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{\chi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{\prime \prime }}\mathit{\delta }{\mathit{v}}^{\mathrm{2}}\mathit{.}$(9)Since we use cubic interpolation for the template, the second derivative of the template F (and of χ²) is continuous. Then a parabolic interpolation through the minimum of the χ²(v) function and the two adjacent neighbours provides a refined estimate for v and an error estimate for v can be obtained from the parabola curvature. The parabola minimum is located at the point, $\mathit{v}\mathrm{=}{\mathit{v}}_{\mathit{m}}\mathrm{-}\frac{\mathrm{\Delta }\mathit{v}}{\mathrm{2}}\mathrm{·}\frac{{\mathit{\chi }}_{\mathit{m}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\chi }}_{\mathit{m}\mathrm{-}\mathrm{1}}^{\mathrm{2}}}{{\mathit{\chi }}_{\mathit{m}\mathrm{-}\mathrm{1}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\chi }}_{\mathit{m}}^{\mathrm{2}}\mathrm{+}{\mathit{\chi }}_{\mathit{m}\mathrm{+}\mathrm{1}}^{\mathrm{2}}}\mathit{,}$(10)where m indexes the discrete global minimum of ${\mathit{\chi }}_{\mathit{k}}^{\mathrm{2}}$ (see Eq. (10.2.1) in Press et al. 1992, that is specialised here for a uniform grid, Δv = v_k−v_k−1, or Eq. (A.3)).

The uncertainty of v is estimated from Eq. (9) with the $\mathrm{\Delta }{\mathit{\chi }}^{\mathrm{2}}\mathrm{(}\mathit{\delta v}\mathrm{=}{\mathit{ϵ}}_{\mathit{v}}\mathrm{)}\mathrm{=}{\mathit{\chi }}^{\mathrm{2}}\mathrm{-}{\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{=}\mathrm{1}$ criterion along with the curvature from Eq. (A.4) ${\mathit{ϵ}}_{\mathit{v}}^{\mathrm{2}}\mathrm{=}\mathrm{2}\frac{\mathrm{1}}{\mathrm{(}{\mathit{\chi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{\prime \prime }}}\mathrm{=}\mathrm{2}\frac{\mathrm{\Delta }{\mathit{v}}^{\mathrm{2}}}{{\mathit{\chi }}_{\mathit{m}\mathrm{-}\mathrm{1}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\chi }}_{\mathit{m}}^{\mathrm{2}}\mathrm{+}{\mathit{\chi }}_{\mathit{m}\mathrm{+}\mathrm{1}}^{\mathrm{2}}}\mathrm{·}$(11)This value might be rescaled with ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ to account for under- or overdispersion of the fit.

Finally, a statistical analysis of the fit is performed. The χ² and ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ are computed. Linear residuals r_i = f_i−f_mod,i exceeding the threshold $\left|{\mathit{r}}_{\mathit{i}}\right|\mathit{>}\mathit{\kappa }{\mathit{ϵ}}_{\mathit{i}}\sqrt{{\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}}\mathit{,}$(12)where the clipping value κ is typically 3...5, are flagged as outliers (e.g. by setting their weights to zero w_i = 0) and excluded in a repeated fitting. The same can be formulated with normalised residuals ${\mathit{\chi }}_{\mathit{i}}\mathrm{=}\frac{{\mathit{r}}_{\mathit{i}}}{{\mathit{ϵ}}_{\mathit{i}}}$$\mathrm{\left|}{\mathit{\chi }}_{\mathit{i}}\mathrm{\right|}\mathit{>}\mathit{\kappa }\sqrt{{\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}}\mathit{.}$(13)Figure 1 illustrates the best fit between one observation and a template.

Each order is fitted separately and finally a weighted mean for the radial velocities v_o from Eq. (10) with errors ϵ_vo from Eq. (11) over all orders o is computed (see also Fig. 3) $\mathit{v}\mathrm{=}\frac{\sum {\mathit{ϵ}}_{{\mathit{v}}_{\mathit{o}}}^{-2}{\mathit{v}}_{\mathit{o}}}{\sum {\mathit{ϵ}}_{{\mathit{v}}_{\mathit{o}}}^{-2}}\mathit{,}$(14)with the error estimate ${\mathit{ϵ}}_{\mathit{v}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\sum {\mathit{ϵ}}_{{\mathit{v}}_{\mathit{o}}}^{-2}}\mathrm{·}\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{o}}\mathrm{-}\mathrm{1}}\begin{array}{c}\sum \end{array}\frac{\mathrm{(}{\mathit{v}}_{\mathit{o}}\mathrm{-}\mathit{v}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{ϵ}}_{{\mathit{v}}_{\mathit{o}}}^{\mathrm{2}}}}\mathrm{=}\frac{\mathrm{wrms}}{\sqrt{{\mathit{N}}_{\mathit{o}}\mathrm{-}\mathrm{1}}}\mathit{,}$(15)where wrms is the weighted root mean square and N_o is the number of orders.

2.3. Spectra coadding

At first glance, coadding spectra seems straightforward. However, the stellar spectra are Doppler shifted by potential Keplerian orbits and by the barycentric motion of Earth. Thus, spectra from different epochs are sampled at different wavelength footpoints. Hence naive coadding would require either some kind of resampling or interpolation of the observations onto a common wavelength grid (e.g. Anglada-Escudé & Butler 2012), which leads to difficulties in applying the data uncertainties, or calculating some form of bin means, which ignores the local gradient over the bin width. We point out that bin means are zeroth-order B-splines.

We carry out the “co-adding” with a uniform cubic basic spline (B-spline) regression to the normalised data. We emphasise here the term regression, which means least-squares fit and should be distinguished from spline interpolation and smoothing splines. There are several benefits of this approach: (1) B-spline regression is a linear least-squares method, and, therefore, it is fast; (2) there is no need to interpolate the data; (3) data point uncertainties can be easily taken into account; (4) robust statistics for outlier detections can be obtained with kappa-sigma clipping; and (5) a spline function is a direct outcome and consistent with our input for the forward model.

Given v_n for the N observations from Sect. 2.2, we recompute the polynomials p_n,i for each order (now with the mean RV v_n) to normalise the data ($\frac{{\mathit{f}}_{\mathit{n,i}}}{{\mathit{p}}_{\mathit{n,i}}}$) and calculate the Doppler-shifted wavelengths with Eq. (5). We factor ${\mathit{p}}_{\mathit{n,i}}^{\mathrm{2}}$ in Eq. (3) and write it in the form ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{n,i}}{\mathit{p}}_{\mathit{n,i}}^{\mathrm{2}}{\mathit{w}}_{\mathit{n,i}}{\left[\frac{{\mathit{f}}_{\mathit{n,i}}}{{\mathit{p}}_{\mathit{n,i}}}\mathrm{-}\mathit{F}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{n,i}}^{\mathrm{\prime }}\mathit{,b}\mathrm{\right)}\right]}^{\mathrm{2}}\mathit{.}$(16)Now the task is to find all the coefficients b_k of the spline with K knots, so that the residuals are minimal. The knots of the template are positioned in a uniform grid in (logarithmic) wavelengths λ_k (lnλ_k) and f_k (or B-spline coefficient b_k, respectively). By default, the number of knots K is similar to the number of data points per spectrum in the order, i.e. we have about one knot per pixel. Depending on the needs, the sampling could also be decreased to obtain a smoother template (e.g. noisy observations or fast rotators) or increased to obtain subpixel sampling (many observations). It can be inferred from Eq. (16) that the normalised error estimates are $\frac{{\mathit{ϵ}}_{\mathit{n,i}}}{{\mathit{p}}_{\mathit{n,i}}}$.

Equation (16) is a linear least-squares problem due to the linearity of the coefficients in the B-spline (de Boor 1978; Dierckx 1993; Eilers & Marx 1996) used to describe the template $\mathit{F}\mathrm{\left(}\mathit{x,b}\mathrm{\right)}\mathrm{=}\begin{array}{c}\sum \end{array}{\mathit{b}}_{\mathit{k}}{\mathit{B}}_{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(17)where B are the basic functions. Figure 2 shows the normalised and Doppler-shifted data and the best spline fit.

However, caution is needed when there are large gaps in the data, i.e. regions where we have no information to say anything about the template (i.e. no constraints for b_k). Possible solutions are (1) splitting the spline fit at those points; (2) fitting penalised splines; or (3) heavily down-weight, but not fully reject, the flagged data points (outliers, tellurics), which caused the gaps. We chose the last option (down-weighting) for tellurics. This means that we include regions heavily contaminated by telluric lines in coadding, while in RV measurements they are totally excluded.

For each knot we estimate an error for the knot value as well as the number of contributing good points. We calculate the errors in the knot values by first estimating the error in the B-spline coefficients as $\begin{array}{ccc}{\mathit{ϵ}}_{{\mathit{b}}_{\mathit{k}}}\mathrm{=}\mathrm{1}/\sqrt{\sum _{\mathit{i}}{\mathit{w}}_{\mathit{i}}{\mathit{B}}_{\mathit{k}}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}}\mathit{,}& & \end{array}$(18)and then through error propagation ${\mathit{ϵ}}_{{\mathit{f}}_{\mathit{k}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{6}}{\mathit{ϵ}}_{{\mathit{b}}_{\mathit{k}\mathrm{-}\mathrm{1}}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{4}}{\mathrm{6}}{\mathit{ϵ}}_{{\mathit{b}}_{\mathit{k}}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{6}}{\mathit{ϵ}}_{{\mathit{b}}_{\mathit{k}\mathrm{+}\mathrm{1}}}^{\mathrm{2}}\mathit{.}$(19)This simplified estimate does not use the covariance matrix, since calculating the inverse matrix would be very time consuming.

To increase the robustness against outliers, a few (≤3) kappa-sigma clipping iterations (κ = 5) are performed for the coadding (similar as in Eq. (13)).

Problems can arise with this coadding method when the pixel phase coverage is small because of a small barycentric range owing to close observations or high ecliptic latitude. Such small coverage may lead to ringing in the spline function. Reverting to simple bin means (zeroth order B-splines) or using penalised splines (p-splines, as a generalisation of B-splines; Eilers & Marx 1996) may help in such situations. Also sharp (undersampled) features, such as cosmics and emission lines, may lead to ringing features. The SERVAL code offers the option to use p-splines, but the choice for the value of penalty/smoothing parameter in each particular case will require some user testing or cross-validation algorithms. We also point out the Gaussian process approach of Czekala et al. (2017), where the variance/smoothness parameter is a hyperparameter, and that B-splines and p-splines can be considered as a special case of Gaussian processes (Rasmussen & Williams 2006).

SERVAL also offers the possibility to use external templates, e.g. from other observed stars of similar spectral type or synthetic spectra. The latter may be useful to derive absolute RVs and line indices (Sect. 3.3).

3. Spectral diagnostic and activity indicators

3.1. Chromatic RV index

Fig. 3 Radial velocities measured in 42 orders of two CARMENES VIS observations of YZ CMi along with the simple weighted average (dashed black line) and the chromatic index (solid blue and red line). The wavelength axis is logarithmic.

Fig. 4 Two Gaussian functions (solid black and dashed red) with slightly different shifts (left panel) and slightly different widths with same amplitude (middle panel) and same unit area (right panel). Their scaled difference (dotted green) is related to the first and/or second derivative (dash-dotted blue).

In Eq. (14) we compute a simple weighted average for the RV over the orders. However, since the echelle orders are related to wavelength, we can also try to get some information about wavelength dependency for instance by simply extending the model to a straight line $\mathit{v}\mathrm{\left(}\mathit{o}\mathrm{\right)}\mathrm{=}\mathit{\alpha }\mathrm{+}\mathit{\beta }\ln {\mathit{\lambda }}_{\mathit{o}}\mathit{,}$(20)where λ_o is a representative wavelength of echelle order o (e.g. lnλ_o = ⟨lnλ_o,i⟩). Via a χ²-fit we obtain a best fit estimate for the slope parameter β, which we call chromatic index. Its unit is velocity per wavelength ratio e (Neper, symbol Np). Figure 3 illustrates the slope definition for two observations of the active M dwarf YZ CMi.

The parameter α is considered as a nuisance parameter. Using the mean velocity v from Eq. (14) and the substitution α = v−βlnλ_v, we can re-parametrise Eq. (20) as $\mathit{v}\mathrm{\left(}\mathit{o}\mathrm{\right)}\mathrm{=}\mathit{v}\mathrm{+}\mathit{\beta }\ln \frac{{\mathit{\lambda }}_{\mathit{o}}}{{\mathit{\lambda }}_{\mathit{v}}}\mathit{,}$(21)where λ_v is the wavelength at which the slope intersects the weighted mean RV. This allows us to report an effective wavelength, which is in particular useful when comparing data from instruments covering different wavelengths.

In the definition of the chromatic index in Eq. (20), we chose lnλ, the natural logarithm of wavelength λ, as independent variable. While other choices like wavelength λ or order o, which is basically equivalent to reciprocal wavelength λ^-1 (i.e. frequency), lead to qualitatively very similar results for β, our choice gives the same weights to each velocity element (resolution element) and is also useful when comparing data from a variety of instruments. Moreover, the wavelength range of 1 Np is typically covered by high resolution echelle spectrographs, i.e. the values of the slopes are of the order of the RV change over all echelle orders.

3.2. Differential line width (dLW)

A feature of the CCF method is that the CCF can be interpreted as a mean stellar line profile that is actually convolved with a kind of kernel. By analysing the CCF shape, for example by fitting a Gaussian function, we can obtain information about the moments, among them the centre of the Gaussian (first moment, i.e. RV) and the width (second moment, i.e. full width half maximum, FWHM). To obtain information about asymmetries (third moment, skewness), often bisectors are used (e.g. Figueira et al. 2013).

We study if there is an analogous set of parameters that can be associated with least-squares fitting. Since the CCF is basically a χ² function, one could consider fitting a Gaussian function to the χ² function. However, we argue that the shape (curvature, asymmetry) of the χ² function is formed by the signal-to-noise of the observation, masking of outliers and tellurics, and also the simultaneous fit of the background polynomial.

In the following, we suggest another approach. We review a simple method to measure differential RVs as implied in Bouchy et al. (2001). Figure 4 (left) shows that when a Gaussian function is slightly displaced, the residuals, i.e. the difference between both curves, are correlated with the first derivative of the Gaussian function. In fact, this is actually the definition of the first derivative. From the scaling factor, which was already applied in Fig. 4, we can derive the RV $\begin{array}{ccc}\mathit{g}\left(\mathit{x}\mathrm{+}\frac{\mathit{v}}{\mathit{c}}\right)\mathrm{-}\mathit{g}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\approx }\frac{\mathit{v}}{\mathit{c}}\hspace{0.17em}{\mathit{g}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}& & \end{array}$Now we replace on the left side the function difference by the residuals r_i between (drifted) data and (flux-scaled) template, and on the right side the derivative of the Gaussian function by the derivative of the template. The derivative must be with respect to velocity, i.e. ${\mathit{f}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\ln \mathit{\lambda }}\mathrm{=}\mathit{\lambda }\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{\lambda }}$. The relation becomes ${\mathit{r}}_{\mathit{i}}\mathrm{=}\frac{\mathit{v}}{\mathit{c}}\hspace{0.17em}{\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(22)In simple terms, we scale the first derivative to the residuals. The best scaling factor v can be derived with a linear least-squares fit weighted with flux uncertainties σ_i. In Sect. 2.2, on the other hand, the RV shift was already applied by actually minimizing this correlation.

Now the question is whether we can find a similar relation using the first or higher derivatives to check for variations in the line width. Therefore, we vary the line width of the Gaussian. In one variant we keep the line height constant (Fig. 4, middle) and in another variant we keep the same area under the curve (Fig. 4, right). For both cases we plot again the scaled residuals and compare these residuals with squared first derivative (middle) and second derivative (right), respectively.

The case of unit area (right panel) results in balanced positive and negative deviations, which is more similar to the situation we are facing after a least-squares fit of the RV (Sect. 2.2). Indeed, we prove in Appendix B that the mean value of the second derivative is zero (the same is true for least-squares fit residuals) and that $\mathit{g}\mathrm{(}\mathit{x,\sigma }\mathrm{+}\mathit{\delta \sigma }\mathrm{)}\mathrm{-}\mathit{g}\mathrm{\left(}\mathit{x,\sigma }\mathrm{\right)}\mathrm{\approx }\mathit{\sigma }\mathrm{·}\mathit{\delta \sigma }\mathrm{·}{\mathit{g}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}$(23)This relation is in a strict sense valid only for Gaussian functions, but it suggests simply scaling the second derivative of the template (${\mathit{f}}^{\mathrm{\prime \prime }}\mathrm{=}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}\ln {\mathit{\lambda }}^{\mathrm{2}}}\mathrm{=}{\mathit{\lambda }}^{\mathrm{2}}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{f}}{\mathrm{d}{\mathit{\lambda }}^{\mathrm{2}}}$) to the residuals ${\mathit{r}}_{\mathit{i}}\mathrm{=}\frac{\mathit{\sigma }\mathrm{·}\mathrm{\Delta }\mathit{\sigma }}{{\mathit{c}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{f}}^{\mathrm{\prime \prime }}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(24)Hence we have to assume that the stellar lines are Gaussian-like shaped. The relation also holds when we build up a spectrum by adding Gaussian lines at other positions and with different strengths due to the linearity. Then this also implies that it is applicable to blended lines. However, the lines should have all the same width. We point out that similar assumptions are made when fitting a Gaussian profile to the CCF.

The scaling factor Δσ carries information about line width changes and we compute its value via a weighted linear least-squares fit resulting in $\mathrm{dLW}\mathrm{\equiv }\mathit{\sigma }\mathrm{\Delta }\mathit{\sigma }\mathrm{=}{\mathit{c}}^{\mathrm{2}}\frac{\sum {\mathit{w}}_{\mathit{i}}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime \prime }}{\mathit{r}}_{\mathit{i}}}{\sum {\mathit{w}}_{\mathit{i}}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime \prime }\mathrm{2}}}\mathit{,}$(25)and estimate its uncertainty through error propagation as ${\mathit{ϵ}}_{\mathit{\sigma }\mathrm{\Delta }\mathit{\sigma }}\mathrm{=}{\mathit{c}}^{\mathrm{2}}\sqrt{\frac{\mathrm{1}}{\sum {\mathit{w}}_{\mathit{i}}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime \prime }\mathrm{2}}}}\mathrm{·}$(26)The dLW is computed in each order and finally averaged similar to the RVs in Eqs. (14) and (15). The SERVAL programme propagates the uncertainties of the spectra which are mostly photon and readout noise. On top of this there can be other noise sources from the instrument, such as focus or resolution change, or observation, such as line broading due to barycentric motion during the exposure. Therefore, even quiet stars with presumable stable intrisic line shapes may show excess variations in line width indicators.

Since our template is a cubic B-spline, we can easily calculate the second derivative f′′ at each position of the spectrum. Also the third derivative would be possible and the possibility of defining a differential alternative to the bisector span in the CCF would appear very straightforward, but this is not yet implemented. Of course, the uncertainty in the higher order moments increase.

Owing to our differential approach, our width indicator σΔσ has also a differential nature, and, therefore, we call it differential line width (dLW). Its unit is m²/ s² if the derivative is calculated in velocity or logarithmic wavelength scale. In analogy to differential RVs, in which an additive offset remains with respect to to absolute RVs, the offset now becomes multiplicative in the differential width (second moment). The σΔσ indicator is sensitive to regions where the second derivative is large, such as the cores of spectral lines (cf. Fig. 4, right panel).

Variations in the spectral line width can be intrinsic to the star, for example pulsation or activity, or of instrumental origin in the form of focus change and smearing due to the barycentric motion during exposure. In any case, correlation of a RV signal with line width variations would argue against a planet hypothesis.

3.3. Line indices

The SERVAL code also provides indices in the form of time series data for a number of spectral lines (e.g. Ca ii H&K, Hα, Na i D, and Ca ii IRT). Before measuring the line indices, we need to find the line positions, i.e. the absolute RVs. These can be obtained by measuring the RV of one spectrum against an absolute reference (e.g. a PHOENIX spectrum).

Following Kürster et al. (2003), the indices are computed as $\mathit{I}\mathrm{=}\frac{\mathrm{⟨}{\mathit{f}}_{\mathrm{0}}\mathrm{⟩}}{\mathrm{0.5}\mathrm{·}\mathrm{(}⟨{\mathit{f}}_{\mathrm{1}}⟩\mathrm{+}⟨{\mathit{f}}_{\mathrm{2}}⟩\mathrm{)}}\mathit{,}$(27)where ⟨f₀⟩ is the mean flux around the line centre and ⟨f₁⟩ and ⟨f₂⟩ are the mean fluxes over reference regions.

As an example, we choose for Hα (λ_air = 6562.8 Å) the region [−40, +40] km s^-1 and as reference [−300, −100] km s^-1 and [+100,+300] km s^-1. We apply a larger core region compared to Kürster et al. (2003), who only applied it to the inactive star GJ 699, to collect all Hα emission for most of the M dwarfs, while for most active stars an even wider range could be considered at the cost of decreasing index precision.

We estimate the error in the mean fluxes using the data error ϵ_i$\mathit{ϵ}⟨\mathit{f}⟩\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}\sqrt{\begin{array}{c}\sum \end{array}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}}\mathit{.}$(28)This choice is not unique, but it preserves some information about the signal-to-noise ratio in the spectrum; in case of pure photon noise ${\mathit{ϵ}}_{\mathit{i}}\mathrm{=}\sqrt{{\mathit{f}}_{\mathit{i}}}$, it satisfies $\mathit{ϵ}⟨\mathit{f}⟩\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{N}}⟨\mathit{f}⟩}$. From the stand point of least squares, the standard error of the mean would be calculated via the standard deviation as $\frac{\mathrm{1}}{\sqrt{\mathit{N}}}\sqrt{\frac{\mathrm{1}}{\mathit{N}\mathrm{-}\mathrm{1}}\sum \mathrm{(}{\mathit{f}}_{\mathit{i}}\mathrm{-}⟨\mathit{f}⟩{\mathrm{)}}^{\mathrm{2}}}$ and for a weighted mean as $\sqrt{\frac{\mathrm{1}}{\sum {\mathit{ϵ}}_{\mathit{i}}^{-2}}}$ (cf. Eq. (15)). However, when modelling a line such as Hα with a mean, i.e. we fit a box, we should be aware of a large model mismatch with this simplistic model making those error estimates misleading. Indeed, one could consider more tailored models, for example the Mount Wilson S-index employs a triangular shape (Duncan et al. 1991).

Finally, we simply propagate the errors from Eq. (28) through Eq. (27) to an uncertainty for the line index ${\mathit{ϵ}}_{\mathit{I}}\mathrm{=}\mathit{I}\sqrt{\frac{{\mathit{ϵ}}_{\mathrm{0}}^{\mathrm{2}}}{{⟨{\mathit{f}}_{\mathrm{0}}⟩}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{2}}^{\mathrm{2}}}{{⟨{\mathit{f}}_{\mathrm{1}}⟩}^{\mathrm{2}}\mathrm{+}{⟨{\mathit{f}}_{\mathrm{2}}⟩}^{\mathrm{2}}}}\mathrm{·}$(29)Sometimes (pseudo-) equivalent widths (pEW) are preferred instead of line indices. We outlined in Zechmeister et al. (2009) that pEWs are closely related to line indices, which might be interpreted as “pseudo line heights”. Indeed, both quantities measure the zero moment (area, integrated flux) and condense it into one parameter. This simple estimate is sufficient for studies on the activity level of stars and temporal variations. However, the true line shape is more complex (e.g. self-absorption features) and needs to be described by more moments/parameters.

We have also implemented a differential version for spectral line indices where the mean fluxes in Eq. (27) are replaced by the scaling factor between the observation and the coadded template. This method aims for optimal weights and highest precision and would be also applicable in case of missing data points, for example due to cosmics or tellurics, but we recommend this only for inactive stars.

4. Results

In this section we evaluate the performance of SERVAL using its RV precision, the dLW, and the chromatic index. To achieve this, we use data from HARPS and CARMENES. The HARPS instrument and its data reduction software has demonstrated long-term 1 m/s precision since its installation in 2003 (Mayor et al. 2003). Each of the points to be tested are evaluated for (1) a stable and inactive M dwarf (Barnard’s star), (2) a G4 dwarf with strong correlations between FWHM and RVs (ζ¹ Ret), and (3) two very active M dwarfs (YZ CMi and GJ 3379). The time-series measurements are shown in Figs. 5−8.

Fig. 5 Radial velocity (left), line width (middle), and chromatic index (right) time series for HARPS data of GJ 699. SERVAL results are indicated in blue solid diamonds and DRS in red open diamonds.

Fig. 6 Same as Fig. 5 but for ζ1 Ret.

Fig. 7 Radial velocity (left), dLW (middle), and chromatic index (right) time series for CARMENES data of YZ CMi.

Fig. 8 Same as Fig. 7, but for GJ 3379 and HARPS (blue diamonds) and CARMENES (green circles) data.

4.1. Targets

GJ699.

Barnard’s star is an inactive M4V star. For comparison we use the same 22 HARPS observations as previously used in Anglada-Escudé & Butler (2012). The spectra were obtained between 2007-04-04 and 2008-05-02. All spectra were secured without simultaneous drift measurements and two spectra ⁴ were taken during twilight.

ζ¹ Ret.

This G4 dwarf was found by Zechmeister et al. (2013) to have a strong correlation between RV and FWHM and hence we selected it as an example for dLW analysis. We analysed 61 HARPS measurements taken during 24 nights over a time span of 1401 days. Typical FWHM values of 7.1 km s^-1 indicate that the star is not a slow rotator.

YZ CMi.

This star (GJ 285, Karmn J07446+035) is a very active M4.5 dwarf with significant rotation measured in several works (vsini ~ 5 km s^-1, log (L_Hα/L_bol) = −3.48; Reiners et al. 2012; Jeffers et al. 2017), a photometric period of 2.78 d (West et al. 2015), an average magnetic field of the order of 3–4 kG, and local field strengths up to 7 kG (Shulyak et al. 2014). We obtained 45 CARMENES VIS observations spanning 480 days (from 2016-01-08 to 2017-05-01).

GJ 3379.

This star (G 099-049, Karmn J06000+027) is an active M4.0 dwarf with significant rotational broadening (vsini = (7.4 ± 0.8) km s^-1, Delfosse et al. 1998). We found 16 HARPS spectra in the ESO archive and obtained 14 spectra with CARMENES for GJ 3379. We only used this star to compare the chromatic index between HARPS and CARMENES wavelength ranges.

4.2. Criterion 1: Radial velocity precision

GJ 699.

Figure 5 (left) compares the RVs derived with SERVAL (rms = 1.30 m/s) and the HARPS DRS pipeline (rms = 1.54 m/s), while Anglada-Escudé & Butler (2012) reported a dispersion of rms = 1.24 m/s for the same dataset. Secular acceleration is subtracted in all cases. This demonstrates the capability of SERVAL to achieve precise RVs at the 1 m/s level.

ζ¹ Ret.

The RVs have an rms of 10 m/s (Fig. 6). Since there is activity-induced jitter, we cannot evaluate the individual RV performance with this star. Still, the RV difference between SERVAL and DRS CCF has an rms of 1.2 m/s and points to 1 m/s performance of SERVAL for solar-like stars as well.

YZ CMi.

The RVs show an rms of 85 m/s (Fig. 7, left). The variation takes place on short time scales (days) with a significant periodicity of 2.776 d with an amplitude of 120 m/s (Fig. 9 and a less significant one-day alias at 1.556 d and 100 m/s). This period also appears in the Hα-line equivalent width but at low significance.

Fig. 9 Radial velocities of YZ CMi phase folded to a period of 2.776 d. The rotational phase is colour coded in reference to Fig. 12; blue and red correspond to the phases of predicted maximum blue shift and red shift, respectively.

4.3. Criterion 2: dLW

GJ 699.

The comparison between dLW and FWHM, drawn on the left and right axes, respectively, is shown as a time series in the middle panel of Fig. 5. Both indicators have relatively similar error bars and vary significantly; we adopted ϵ_FWHM = 2.35ϵ_RV. Using the relation FWHM = 2.35σ for a Gaussian, we would expect $\mathrm{\Delta }\mathit{FWHM}\mathrm{=}\frac{\mathit{FWHM}\mathrm{·}\mathrm{\Delta }\mathit{FWHM}}{\mathit{FWHM}}\mathrm{=}\frac{{\mathrm{2.35}}^{\mathrm{2}}}{\mathit{FWHM}}\mathrm{·}\mathit{\sigma }\mathrm{\Delta }\mathit{\sigma }$. This theoretical prediction is shown as a blue line in Fig. 10 (top left panel). We see that dLW values scatter around this prediction, but the large uncertainty for the best fitting slope shows that the dLW and FWHM correlate poorly for GJ 699. Surprisingly, we find a much better correlation of dLW with FWHM divided by the contrast indicator, which is the relative amplitude of the Gaussian fit (Fig. 10, top right). This is currently not well understood, but it likely reminds us that dLW and FWHM are not fully tracking the same effect. Indeed, the dLW indicator assumes that the area of the Gaussian-shaped line, i.e. the product of FWHM and contrast, remains constant, while a Gaussian fit to the CCF decomposes both parameters simultaneously. Hence, we could argue that additional contrast variations (on top of contrast ∝ FWHM^-1) could lead to biased dLW measurements. Those contrast variations could be induced by imperfect background subtraction or have stellar origin, albeit we believe the latter is less likely for this quiet star. Alternatively, the poor dLW-FWHM correlation for GJ699 could point to a bias when fitting a Gaussian function to the CCF of an M dwarf, which is known to exhibit prominent side lobes (see e.g. Fig. 7 in Berdiñas et al. 2017). Neither dLW nor FWHM correlate with RVs.

Fig. 10 Left: correlation between the DRS-FWHM and the SERVAL-dLW with HARPS data for the stars GJ 699 (top) and ζ1 Ret (bottom). The solid black line indicates the best fitting linear trend; the blue dashed line represents the theoretical prediction. Right: same as left panels, but the FWHM is divided by the contrast.

ζ¹ Ret.

The time series of dLW and FWHM already indicates a similar behaviour (Fig. 6, middle). Both, FWHM and dLW, correlate with RV (Fig. 11). There is also a direct correlation between dLW and FWHM; still the correlation improves again when dividing the FWHM by the contrast (Fig. 10, bottom panels).

Fig. 11 Correlation between RV and dLW for HARPS data of ζ1 Ret.

YZ CMi.

The dLW (Fig. 7, middle) reveals a 2.776 d period that we have already found in the RVs. The correlation between dLW and RV is not linear (Fig. 12, bottom and top), but we find that our observations follow a well-defined path when we colour code the rotational phase with this period. The one dLW outlier is an observation with an Hα flare. The spectrum has the highest Hα core emission and noticable rising Hα wings compared to the other spectra. Likely, the flare adds continuum flux over the full spectrum, leading primarily to line contrast variations. The flare event is not noticable in RV and chromatic index, although the point has maximum RV and minimum chromatic index.

4.4. Criterion 3: chromatic RV Index

GJ 699.

The chromatic index (Fig. 5, right) does not vary significantly. This is expected for an inactive and slowly rotating star.

ζ¹ Ret.

The chromatic index (Fig. 6) shows long-term variations that seem to be anti-correlated with the RV long-term variations, which could be an activity cycle. However, the chromatic index does not follow the short-term variations of RV and dLW. Therefore, the chromatic index apparently has a complementary role to line width indicators.

Fig. 12 Top: chromaticity (correlation between RV and chromatic index) for 45 CARMENES VIS spectra of YZ CMi. The left-most (RV = 165.3 m/s, β = 453 m/s/Np) and right-most data points (RV = 136.6 m/s, β = −441 m/s/Np) are the values derived in Fig. 3. Bottom: RV-dLW correlation is shown.

YZ CMi.

In the top panel of Fig. 12, we show the chromatic index as a function of RV for our CARMENES observations of YZ CMi; both parameters are clearly anti-correlated. We call this correlation and the corresponding slope chromaticity κ. RV variations of ± 150 m/s correspond to variations of the chromatic index of ∓ 450 m/s/Np. Of course, due to the anti-correlation, the periodograms for chromatic index and RV have a similar shape and the 2.776 d period is also found in the chromatic index. For this period, we colour code our data in Fig. 12 according to their rotational phase. One can clearly follow how the individual observations are distributed throughout the rotational phase of YZ CMi.

The chromatic index-RV correlation demonstrates that the dominating effect in the RV variations of YZ CMi is wavelength dependent, which is consistent with the hypothesis that they are caused by co-rotating features such as active regions on the surface of the star. The wavelength-dependent temperature contrast between quiet and active regions is expected to cause a variation of the RVs as a function of wavelength (e.g. Fig. 12 in Reiners et al. 2010). Active regions (or spots) that are somewhat cooler (or hotter) than the photosphere show less contrast to the photosphere at longer wavelengths. The chromatic index attempts to capture this dependence on wavelength. The advantage of the chromatic index over chromospheric emission (Ca ii H&K, Hα) is that it comes with a directional sense, i.e. it can be positive and negative. The chromatic index is directly connected to RV because it is calculated from the same line profile deformation.

A negative chromaticity, as in case of YZ CMi, means that the RV scatter decreases towards redder wavelength. This amplitude decrease can been seen in Fig. 3 and is predicted with spot simulations (Desort et al. 2007; Reiners et al. 2010). The slope in the correlation plot (Fig. 12, top) provides information about spot temperatures. Moreover, we identify a separation between the values of the chromatic index at rotational phases 0.5 and 1, which would not be the case if the relation was simply linear and is probably related to convective blue shift. A deeper interpretation of these effects requires more detailed observations and goes beyond the scope of this paper.

Using the example of YZ CMi, we demonstrated that the chromatic index in the wavelength range of CARMENES VIS (550–970 nm) is a powerful tool for understanding RV variations induced by stellar activity. A closer look at the individual RVs per wavelength (Fig. 3) shows that a fairly tight relation between RV and wavelength exists for the wavelength range 600–920 nm, but the five spectral orders short of λ = 590 nm fail to follow this trend. We observe this behaviour in the other exposures of YZ CMi, too. The significance of this effect for YZ CMi and other stars needs more detailed investigation. There are two plausible explanations for the lack of correlation at shorter wavelengths: first, the contrast grows too high in that active regions no longer contribute to the observed spectrum in a significant way; and second that the intensity of spectral features at short wavelengths also differs dramatically between active and quiet regions (Reiners et al. 2010).

A possible way to investigate whether the chromatic index effect continues down to shorter wavelengths is to determine RVs for individual spectral orders of spectra taken with HARPS. This was attempted for the active M4.5 dwarf AD Leo by Reiners et al. (2013, Fig. 9). Interestingly, the activity-induced RV amplitude reported there grows with wavelength (positive chromaticity κ), but the effect is more than an order of magnitude weaker than what we observe in the CARMENES data of YZ CMi. Additionally, the two RV amplitudes shown for AD Leo by Reiners et al. (2013) that overlap with the CARMENES wavelength range in fact show a much steeper wavelength dependence with the same negative sign as observed in the CARMENES data for YZ CMi. The comparison between AD Leo and YZ CMi, however, remains inconclusive because we might just be seeing very different effects in different stars.

In order to make a more useful comparison of the HARPS and CARMENES wavelength ranges, we looked for targets that were observed with both instruments, albeit not simultaneously. One of the few targets available for this exercise is GJ 3379. Radial velocities and line indicators are shown in Fig. 8. The datasets were taken at very different times. The RVs from HARPS spectra exhibit an rms of 74 m/s while the CARMENES RVs show an rms of 21 m/s. The intrinsic uncertainties of the individual RV measurements are of the order of a few m/s in both cases and negligible for this comparison.

We calculated the chromatic index in both datasets and show the correlations between chromatic index and RV in Fig. 13. For our CARMENES data, chromatic index and RVs are correlated. The relation can be described by a linear fit with a gradient of −2.6 Np^-1(± 21%). In the HARPS data, we find no correlation between chromatic index and RV, although the scatter in RV is higher than in the CARMENES data.

This comparison seems to indicate that the RVs of GJ 3379 show larger scatter at bluer wavelengths and that the RV variations at longer wavelengths are correlated with the chromatic index, while those at bluer wavelengths are not. However, an alternative explanation is that the RV jitter of the star changed during the course of a potential activity cycle; the HARPS and CARMENES data are separated in time by four years. One way to test this scenario is to compare RVs from the same wavelength range in both datasets; HARPS and CARMENES overlap in the range 550–690 nm. Unfortunately, this wavelength range is too short to draw definite conclusions about the chromatic index, and we have already shown above that the slope is most useful at longer wavelengths. Nevertheless, we find a correlation slope between chromatic index and RV of about −1 Np^-1. Although the scatter around this relation is relatively high, HARPS and CARMENES data points are consistent with each other. The correlation between chromatic index and RVs is much weaker than for the full CARMENES wavelength range and the correlation itself is marginally significant. On the other hand, the RVs calculated from the limited wavelength range are precise enough to compare the two time series; in this case we find that the rms of CARMENES RVs (30 m/s) is significantly lower than the HARPS rms (77 m/s). It is unlikely that this effect is caused by differences between the two instruments because for both we estimate the internal uncertainties to be much smaller and comparable with each other (4 m/s). Thus, for now the case remains undecided. We find that GJ 3379 shows reduced activity jitter during the epoch observed with CARMENES, while the RV rms was high when it was observed with HARPS. Nevertheless, the CARMENES RVs show a correlation with chromatic index when the full CARMENES wavelength range is taken into account. We expect that such a correlation would still be visible when GJ 3379 exhibits larger RV variations, but simultaneous observations of the HARPS and CARMENES wavelength ranges are required to settle this question.

Fig. 13 Radial velocity-chromatic index correlation for GJ 3379. The fitted slope is the chromaticity of GJ 3379 in HARPS data (blue diamonds) and CARMENES data (green circles) using their full spectral range (left) and only the overlapping spectral region (right).

5. Summary

We have presented our concept to derive high-precision RVs and have provided details on the creation of a template. Following an iterative and sequential approach we first derive approximate RVs measured against an observed spectrum and then improve the template by co-adding all observed spectra and recompute the RVs. The spectra are “coadded” via a weighted least-squares regression with a cubic B-spline.

The SERVAL code yields an output consisting of time series for high-precision RVs and a number of spectral indicators useful for further diagnostics as well as the high signal-to-noise templates partly cleaned from telluric contamination taking advantage of barycentric shifts. We find a similar performance when comparing the RV results from SERVAL and the HARPS-DRS pipeline using data for GJ 699 and ζ¹ Ret. An additional example can also be found in Zechmeister et al. (2014), where an older version of SERVAL was applied to ~150 τ Cet spectra obtained during one night. This demonstrates that SERVAL produces RVs at a 1 m/s precision level.

Furthermore, we motivated a definition to study differential changes in the spectral line widths. We scale the second derivative of the template to the residuals from the fit for the best RVs and have shown that it can be useful in practice. This differential indicator nicely fits in our self-consistent framework of least-squares fitting (for differential RVs) without the need for external templates. Of course, this comes at the price that the differential results have no external accuracy, i.e. the line width remains unknown. In contrast, unbiased and accurate line width measurements for M dwarf spectra with their overall blended lines can be obtained with least-squares deconvolution (Barnes et al. 2012), but those are more complicated algorithms that require accurate line positions as input.

The chromatic index, while very simple in its definition, turns out to be a very powerful indicator to identify stellar activity in RV signals. We have shown that it correlates very well with RVs for the active M dwarf YZ CMi. Hence it will help us to identify signals that are associated with the intrinsic stellar variations rather than variations induced by Keplerian motion, to lower the detection limit for active and inactive stars because it can be used to remove and disentangle the activity contribution from the RVs, and to better understand the physics of the stars, such as spot temperature or convective blue shift. The chromatic index summarises in one value a first-order effect of wavelength dependence in the RVs and is well suited for an efficient analysis. Still, a more detailed investigation in smaller wavelength ranges might further our understanding of the wavelength dependence of stellar radial velocities (Reiners et al. 2013; Feng et al. 2017). A large wavelength range, such as that offered by CARMENES, is essential for this.

Acknowledgments

We thank J. Zhao, Y. Thiele, E. Nagel, L. Nortmann, and G. Anglada-Escudé for software testing and helpful discussions. M.Z. acknowledges support from the Deutsche Forschungsgemeinschaft under DFG RE 1664/12-1 and Research Unit FOR2544, project No. RE 1664/14-1. I.R. acknowledges support by the Spanish Ministry of Economy and Competitiveness (MINECO) and the Fondo Europeo de Desarrollo Regional (FEDER) through grant ESP2016-80435-C2-1-R, as well as the support of the Generalitat de Catalunya/CERCA programme. V. Wolthoff acknowledges funding from the DFG Research Unit FOR2544, project No. RE 2694/4-1. V.J.S.B. are supported by grant AYA2015-69350-C3-2-P from the Spanish Ministry of Economy and Competiveness (MINECO). The UCM, CAB, and IAA team members acknowledges support by the Spanish Ministry of Economy and Competitiveness (MINECO) from projects AYA2016-79425- C3-1,2,3-P. CARMENES is an instrument for the Centro Astronómico Hispano-Alemán de Calar Alto

(CAHA, Almería, Spain). CARMENES is funded by the German Max-Planck-Gesellschaft (MPG), the Spanish Consejo Superior de Investigaciones Científicas (CSIC), the European Union through FEDER/ERF funds, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Insitut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the Spanish Ministry of Economy, the German Science Foundation (DFG) through the Major Research Instrumentation Programme and DFG Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía.

References

Appendix A: Parabolic interpolation

We want to derive the interpolating parabola y = ax² + bx + c through the three points (x_i,y_i), (x_i−1,y_i−1), and (x_{i + 1},y_{i + 1}). The points should be equidistant, i.e. x_{i + 1} − x_i = Δx and x_i − x_i−1 = Δx. Now we centre points to (x_i,y_i) and the three transformed points are (0,0), (− Δx,Y_-1) and $\mathrm{\left(}\mathrm{\Delta }\mathit{x,Y}{\mathrm{+}\mathrm{1}}_{}\mathrm{\right)}$. Inserting into the parabola equation leads to a system of three equations $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \mathit{c}\\ {\mathit{Y}}_{-1}& \mathrm{=}& \mathit{a}\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}\mathrm{-}\mathit{b}\mathrm{\Delta }\mathit{x}\\ {\mathit{Y}}_{\mathrm{1}}& \mathrm{=}& \mathit{a}\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}\mathrm{+}\mathit{b}\mathrm{\Delta }\mathit{x}\mathit{.}\end{array}$(A.1)The solution is $\begin{array}{ccc}\mathit{a}& \mathrm{=}& \frac{{\mathit{Y}}_{-1}\mathrm{+}{\mathit{Y}}_{\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}}\mathrm{=}\frac{{\mathit{y}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{-}\mathrm{2}{\mathit{y}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}}\\ \mathit{b}& \mathrm{=}& \frac{{\mathit{Y}}_{\mathrm{1}}\mathrm{-}{\mathit{Y}}_{-1}}{\mathrm{2}\mathrm{\Delta }\mathit{x}}\mathrm{=}\frac{{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{y}}_{\mathit{i}\mathrm{-}\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }\mathit{x}}\\ \mathit{c}& \mathrm{=}& \mathrm{0.}\end{array}$(A.2)Reforming the parabola equation to $\mathit{y}\mathrm{=}{\mathit{y}}_{\mathit{i}}\mathrm{+}\mathit{a}\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathit{i}}\mathrm{+}\frac{\mathit{b}}{\mathrm{2}\mathit{a}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}\frac{\mathit{b}}{\mathrm{2}\mathit{a}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathit{c}$, the extremum of the parabola is at ${\mathit{x}}_{\mathit{c}}\mathrm{=}{\mathit{x}}_{\mathit{i}}\mathrm{-}\frac{\mathit{b}}{\mathrm{2}\mathit{a}}\mathrm{=}{\mathit{x}}_{\mathit{i}}\mathrm{-}\frac{\mathrm{\Delta }\mathit{x}}{\mathrm{2}}\mathrm{·}\frac{{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{y}}_{\mathit{i}\mathrm{-}\mathrm{1}}}{{\mathit{y}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{-}\mathrm{2}{\mathit{y}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}}\mathit{,}$(A.3)and the second derivative at x_c is ${\mathit{y}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{2}\mathit{a}\mathrm{=}\frac{{\mathit{y}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{-}\mathrm{2}{\mathit{y}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}}\mathit{,}$(A.4)which is the definition of the finite second derivative.

Appendix B: Derivatives of the Gaussian function

We start with a Gaussian function with unit amplitude. Their first and second derivatives (with respect to x) are given by $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathit{,}& & \end{array}$(B.1)$\begin{array}{ccc}{\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathit{x}}{{\mathit{\sigma }}^{\mathrm{2}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathrm{=}\mathrm{-}\frac{\mathit{x}}{{\mathit{\sigma }}^{\mathrm{2}}}\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}\\ {\mathit{f}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}}\left(\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}\end{array}$We seek to understand how the residuals should appear when we change the width while keeping the amplitude. The derivative of Eq. (B.1) with respect to parameter σ can be written as $\begin{array}{ccc}\frac{\mathit{\partial f}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{\partial \sigma }}& \mathrm{=}& \frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{3}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathrm{=}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{3}}}\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ & \mathrm{=}& \mathit{\sigma }\frac{{\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}{\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{f}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{·}\end{array}$(B.4)When we are interested in the residuals while conserving the area under the curves (Fig. 4, right panel), we have to include the normalising factor $\frac{\mathrm{1}}{\mathit{\sigma }}$$\begin{array}{ccc}\mathit{g}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{\sigma }}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(B.5)The relation of g′(x) and g′′(x) to g(x) is the same as for f(x). However, the derivative with respect to σ is $\begin{array}{ccc}\frac{\mathit{\partial g}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{\partial \sigma }}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\sigma }}\left(\frac{{\mathit{x}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right)\mathit{g}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ & \mathrm{=}& \end{array}$(B.6)This suggests that we correlate the second derivative of the spectrum with the residuals when we want to infer differential variations in the width.

The mean value of the second derivative is ${\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathit{g}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}\mathrm{=}{\mathit{g}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}|\begin{array}{c}\mathrm{\infty }\\ \mathrm{-}\mathrm{\infty }\end{array}\mathrm{=}\mathrm{0.}$(B.7)

All Figures

	Fig. 1 Least-squares RV fit for RV with CARMENES VIS data for GJ 1265. Top: one observation (red) and the best fit model (blue) are shown. Data points within telluric regions (grey area) are a priori excluded. Data points classified as outliers are indicated in cyan. Bottom: fit residuals are shown.
In the text

	Fig. 2 Spectra coadding of seven CARMENES VIS spectra of GJ 1265. The seven observations (colour coded), each normalised by a background polynomial, are fitted by a uniform cubic B-spline (black) where basically each knot value is a free parameter (black dots). Points affected by tellurics or outlying by 5σ are excluded from the fit (open circles).
In the text

	Fig. 3 Radial velocities measured in 42 orders of two CARMENES VIS observations of YZ CMi along with the simple weighted average (dashed black line) and the chromatic index (solid blue and red line). The wavelength axis is logarithmic.
In the text

	Fig. 4 Two Gaussian functions (solid black and dashed red) with slightly different shifts (left panel) and slightly different widths with same amplitude (middle panel) and same unit area (right panel). Their scaled difference (dotted green) is related to the first and/or second derivative (dash-dotted blue).
In the text

	Fig. 5 Radial velocity (left), line width (middle), and chromatic index (right) time series for HARPS data of GJ 699. SERVAL results are indicated in blue solid diamonds and DRS in red open diamonds.
In the text

	Fig. 6 Same as Fig. 5 but for ζ1 Ret.
In the text

	Fig. 7 Radial velocity (left), dLW (middle), and chromatic index (right) time series for CARMENES data of YZ CMi.
In the text

	Fig. 8 Same as Fig. 7, but for GJ 3379 and HARPS (blue diamonds) and CARMENES (green circles) data.
In the text

	Fig. 9 Radial velocities of YZ CMi phase folded to a period of 2.776 d. The rotational phase is colour coded in reference to Fig. 12; blue and red correspond to the phases of predicted maximum blue shift and red shift, respectively.
In the text

	Fig. 10 Left: correlation between the DRS-FWHM and the SERVAL-dLW with HARPS data for the stars GJ 699 (top) and ζ1 Ret (bottom). The solid black line indicates the best fitting linear trend; the blue dashed line represents the theoretical prediction. Right: same as left panels, but the FWHM is divided by the contrast.
In the text

	Fig. 11 Correlation between RV and dLW for HARPS data of ζ1 Ret.
In the text

	Fig. 12 Top: chromaticity (correlation between RV and chromatic index) for 45 CARMENES VIS spectra of YZ CMi. The left-most (RV = 165.3 m/s, β = 453 m/s/Np) and right-most data points (RV = 136.6 m/s, β = −441 m/s/Np) are the values derived in Fig. 3. Bottom: RV-dLW correlation is shown.
In the text

	Fig. 13 Radial velocity-chromatic index correlation for GJ 3379. The fitted slope is the chromaticity of GJ 3379 in HARPS data (blue diamonds) and CARMENES data (green circles) using their full spectral range (left) and only the overlapping spectral region (right).
In the text