© ESO 2020

1. Introduction

The Doppler spectroscopy or radial velocity (RV) technique is one of the main methods to detect, confirm, and characterise planetary companions. It pioneered the discovery of the first exoplanet orbiting a main sequence star (Mayor & Queloz 1995). Since then, Doppler spectroscopy and photometric searches for transits have become the most successful techniques to detect exoplanets, with about 900 and 3000 discoveries, respectively, of the more than 4000 exoplanets confirmed to date¹.

Early RV searches for exoplanets mainly focused on solar-type stars and, due to their intrinsic faintness, only included a handful of M-dwarf stars. Out of approximately the first hundred exoplanet detections, only a few were found orbiting M dwarfs (Marcy et al. 1998, 2001; Delfosse et al. 1998; Butler et al. 2004). Nevertheless, in the last two decades, interest has rapidly grown regarding M dwarfs as exoplanets hosts, which now account for over 200 exoplanets orbiting about 100 hosts (Martínez-Rodríguez et al. 2019). Compared to more massive and hotter FGK-type stars, the star-to-planet mass ratio is lower for M dwarfs and their habitable zones (HZs) are located closer in. Therefore, the gravitational pull of an HZ-orbiting planet is stronger for M dwarfs, implying RV signals with larger amplitudes. Also, orbital periods are shorter, thus reducing the monitoring timescales. These circumstances make the detection of small, rocky exoplanets feasible with the current instrumentation. Moreover, M dwarfs are the most abundant stars in the solar neighbourhood, which makes them the most common potential exoplanet hosts. The main issue with M dwarfs is that they can exhibit high levels of stellar activity. This activity manifests itself in the form of features in the stellar surface that create distortions on the spectral line profiles used to measure RVs, which complicates the detection of exoplanets. They are also relatively fainter than more massive stars and emit most of their bulk energy in the near-infrared.

The traditional method of measuring Doppler shifts in spectra is to directly cross-correlate them with a template. This is known as the cross-correlation function (CCF) method (Queloz 1995; Baranne et al. 1996). Binary masks are commonly used templates and they represent a very simple model of the stellar spectrum. They consist of a set of boxcar functions (a rectangle-shaped function that has a constant value over a specific range and is zero elsewhere) centred at the minimum wavelength of a large number of absorption lines that are present in the stellar spectrum. In fact, early echelle spectrographs used physical masks, that is, spatial filters with real holes at the positions of absorption lines (Baranne et al. 1979), thus explaining the name. Depending on how much RV information (Bouchy et al. 2001) each line contains, the boxcar functions can have different weights, which improve the RV measurements (Pepe et al. 2002). Deep and narrow lines, that is, lines with the steepest profiles, have more information than shallow and broad lines, hence their weight in the mask is larger.

The CCF method has been the standard approach to obtain precise RV measurements of FGK-type stars. Their spectra display a number of resolved and unblended absorption lines across the visible wavelength range, which are ideal features to be represented by a weighted binary mask. M dwarfs also show a large number of absorption lines. However, due to their lower temperatures, the number of lines is so high that most of the features are blended, making their identification and, hence, the construction of a mask, more challenging. Different methods consisting in a least-squares matching of a full spectrum template to the observations (Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015; Zechmeister et al. 2018), rather than only using a selected set of absorption lines, have proven to work better for M-dwarf stars since they use more Doppler information contained in the spectrum. Recently, different methods have been proposed with the goal of improving RV extraction, such as using individual absorption lines (Dumusque 2018) or approaches using Gaussian processes (Rajpaul et al. 2020).

The CCF continues to be useful for studying the stellar activity of stars, including M dwarfs. The CCF profile gathers information from each individual spectral line selected in the mask, which acts as a kernel, and converts this into an average high signal-to-noise ratio (S/N) line profile. Therefore, average changes in the individual line profiles due to stellar activity are reflected in changes in the CCF shape. Several parametrisations of the CCF profile describing its width, depth, or asymmetry are commonly used proxies of stellar activity. Another benefit is that one can estimate the absolute RV value of the observed stars using RV measurements obtained from the CCF. CCF masks are usually calibrated to an absolute RV scale, while template matching techniques provide RVs that are relative to the template used, which is not necessarily on an absolute scale. Of course, one could calibrate the template against an absolute reference, but the most straightforward way to do that would be by cross-correlation.

In this article, we study the CCF and its activity indicators in M-dwarf stars in the visible and near-infrared spectral range. We use spectral data obtained with the CARMENES instrument (Calar Alto high-Resolution search for M dwarfs with Exo-earths with Near-infrared and optical Echelle Spectrographs, Quirrenbach et al. 2016, 2018), located at the 3.5 m telescope at Calar Alto Observatory in Almería, Spain. It consists of a pair of cross-dispersed, fibre-fed echelle spectrographs with different wavelength coverage: The visual (VIS) channel covers the spectral range λ = 520−960 nm at a resolution of R = 94 600, and the near-infrared (NIR) channel covers the range λ = 960−1710 nm at a resolution of R = 80 400. The CARMENES survey (also guaranteed-time observations – GTO program) has been ongoing since 2016. It monitors over 300 M dwarfs across all spectral subtypes with the main goal of detecting orbiting exoplanets (Reiners et al. 2018a).

In Sect. 2 we present an overview of the CCF method and its parameters. Then we explain the methodology we follow to build binary masks in Sect. 3 and to compute CCFs and their main parameters in detail in Sect. 4. Next, in Sect. 5 we apply these methods to the M dwarf observations of the CARMENES GTO sample. There, we also comment on the masks and CCFs obtained for different targets, compare the RVs obtained with the default CARMENES template matching pipeline, show the behaviour of the CCF activity indicators for an active star in the sample, YZ CMi, and estimate the absolute RV of the CARMENES sample stars using the values derived from their CCFs. Finally in Sect. 6 we summarise our work.

2. Cross-correlation with a weighted binary mask

Measuring the radial velocity of an observed spectrum via cross-correlation relative to a binary mask implies Doppler-shifting the mask to cover a range of velocities when in search for the best match². At each velocity step, the spectrum is multiplied by the mask. Since the mask has a zero value everywhere except at the positions of the selected absorption lines, only the spectrum pixels that totally or partially overlap with a mask line (non-zero values) contribute to the CCF value. For a given velocity shift v, the cross-correlation of a spectrum consisting of n pixels with a mask with m lines is computed according to the following equation

$$\begin{array}{c}\hfill \mathrm{CCF}(v)={\displaystyle \sum _{l=1}^m\sum _{x=1}^n{w}_l{f}_x{\mathrm{\Delta }}_{xl}(v),}\end{array}$$(1)

where f_x is the flux of the spectrum at pixel x, w_l is the weight of the mask line l, and Δ_x l is the fraction of pixel x covered by the mask line l when the mask is shifted by v. This equation corresponds to the cross-correlation function (CCF) of two discrete functions. Their product is computed pixel by pixel, and the term Δ_x l is employed to take partial overlaps between the pixels of the spectrum and the mask into account. The mask lines are defined as having a width of 1 pixel, which corresponds to the minimum interval containing relevant information. The minimum of the CCF corresponds to the Doppler shift of the stellar spectrum with respect to the mask. An analytic function, most frequently a Gaussian, is fitted to the CCF profile for a robust estimate of the position of the minimum and to characterise its profile.

Compared to template matching algorithms, the CCF method can be considered to be a simplistic least-squares minimisation in which we minimise the product of the observed spectrum and model (i.e. the CCF mask); afterwards, we fit the CCF with a function, such as a Gaussian, to estimate the Doppler shift. This is equivalent to computing the CCF between the spectrum and the mask previously convolved with a Gaussian (Baluev & Shaidulin 2015). The only parameter in the model is the RV shift applied to the mask. Other corrections, such as normalisations or slope changes in the spectrum continuum, which are usually forward-modelled in template matching procedures, are not simultaneously modelled when computing the CCF. But rather, they are adjusted a priori. Despite the number of simplifications, the CCF method is capable of yielding precise results. Another difference is that in template matching schemes that use co-added templates, the noise of each of the individual spectra used is present in the template. This creates self-biases towards zero velocity due to the noise correlating with itself, specially if the S/N of the template is low (van Kerkwijk et al. 1995), which should not affect the CCF method.

As mentioned above, the CCF represents an average of the profiles of all the stellar absorption lines that are present in the mask used in velocity space. The different processes in the stellar atmosphere that affect the average position and shape of the individual spectral lines are reflected in the CCF. Therefore, the shape of the CCF contains information on the line profile distortions caused, for example, by stellar activity features present in the photosphere, such as dark spots or faculae.

Parameters, such as the full-width-at-half-maximum (FWHM) and the contrast of the CCF, which are obtained from the Gaussian fit are used to measure temporal changes. Correlations of these parameters with the RV, or the presence of the same periodic signal in the RV and these parameters, may indicate the presence of activity-induced variations in the RVs. They can also trace instrumental variations, such as focus changes, or observational effects, such as barycentric broadening. The bisector of the CCF is used to analyse its asymmetry, which is also an indicator of changes caused by stellar activity for lines of a different strength, correlated with the atmospheric depth of their formation (e.g. Gray 2009). Different ways to quantify the bisector have been proposed. A commonly used metric is the bisector inverse slope (BIS), which was introduced by Queloz et al. (2001). It is defined as the difference between the average velocity of the top region of the CCF (e.g. from 60 to 90%) and the average velocity of the bottom region (e.g. from 10 to 40%). Different fractions of the CCF can be used to compute the BIS, such as the BIS⁺, which uses the 80−90% and 10−20% regions, and the BIS⁻, which uses the 60−70% and 30−40% regions as defined in Figueira et al. (2013). These alternative definitions have been shown to provide more significant results than the usual regions for some specific cases.

Other ways to quantify the asymmetry have also been proposed. In Boisse et al. (2011), the asymmetry was measured as the difference between the RVs obtained by fitting a Gaussian to the upper and lower part of the CCF. This indicator, known as velocity span, V_span, seems to be more robust than the BIS at a low S/N. Figueira et al. (2013) also defined two different indicators. Instead of fitting a Gaussian function to the CCF, they fitted a bi-Gaussian function (two half-Gaussians with different FWHM, each modelling one side of the CCF). This fit of a bi-Gaussian gives an estimate of the RV (the minimum of the CCF) that is different than the usual Gaussian fit, so the difference between these two RVs, ΔV (first used in Nardetto et al. 2006), can be used as an asymmetry measure. In some cases, the amplitude of this indicator is larger than that of the BIS, allowing for the detection of correlations with the RV with a smaller amplitude. The second indicator, V_asy, which was later modified by Lanza et al. (2018), compares the RV information content (the total flux and slope) between each side of the CCF to measure its asymmetry, and it is also found to provide stronger correlations with the RV than other indicators. Simola et al. (2019) also fitted a non-symmetric function to the CCF, in this case a skew normal distribution (a Gaussian with a skewness parameter), to account for asymmetry in the CCF profile. The parameters derived from this fit seem to be more sensitive to stellar activity. Even though different asymmetry parametrisations are proposed in the literature, here, we focus on the original BIS definition because it is the most commonly used asymmetry metric (e.g. it is provided by the HARPS/HARPS-N data reduction software) and we do not find significant differences with the other definitions presented. A more in-depth study of the different asymmetry parametrisations is out of the scope of this work.

3. Weighted binary mask creation method

To compute a CCF, we need to obtain the appropriate weighted binary mask. The main issue is to select a set of absorption lines in the spectra that are appropriate for unbiased radial velocity determination. The lines that we find in the stellar spectrum both mainly depend on the instrument used (i.e. its wavelength coverage and resolution) and the target properties (i.e. its effective temperature and projected rotational velocity). In the case of M dwarfs, the spectrum contains a large number of lines and is rich in molecular bands. The lines are often not clearly separated and, in many cases, they can be severely blended.

Since we want to compute CCFs of M dwarfs observed with the CARMENES instrument, we chose to directly use the large number of available observations from the GTO survey to identify and select lines (Reiners et al. 2018a). We could have also used a high S/N, high-resolution synthetic stellar templates, or information from line lists, but initial tests showed that using our observations yielded clearly better results and the procedure is more straightforward. The number of observations available for all the CARMENES survey targets is sufficient to build high S/N templates by co-adding them.

In the sections below, we describe the procedure we followed to build a mask that is suitable to compute the CCFs. The methods to build masks and compute CCFs with their moments described in the current and next sections have been implemented in the code RACCOON (Radial velocities and Activity indicators from Cross-COrrelatiON with masks). Currently, the code works for CARMENES and HARPS/HARPS-N (Mayor et al. 2003; Cosentino et al. 2012) data, and it is publicly available³.

3.1. High S/N stellar spectrum template

We employed the CARMENES survey observations as the starting point to select the lines. To minimise the effect of the photon noise of individual spectra, we created a high S/N spectral template by co-adding all of the different observations available of a specific star. The templates were obtained using the SERVAL software (Zechmeister et al. 2018), which is the main pipeline used to analyse RVs in the CARMENES survey. SERVAL creates a high S/N stellar template by co-adding observations and then it uses this template to compute a least-squares fit to each of them from which the RV time series is obtained iteratively. Before co-adding, the observations are corrected for the corresponding barycentric motion of the Earth and any other known drift so that the stellar lines are optimally aligned. The templates have a similar format as the observations, and echelle orders are considered individually. The CARMENES NIR orders have a middle gap because of the NIR detector configuration (it is a mosaic of two detectors with a middle gap). Due to this, throughout all this work, we treated the NIR half-orders as single ones.

We then normalised the template orders by dividing each of them by a fit to their “continuum”. To measure the continuum of each order, we selected the pixels with the maximum flux over windows of a few tens to a few hundreds of pixels in width to avoid overfitting as there are numerous absorption lines present in the spectrum. When selecting these pixels, we avoided spectral regions containing very strong lines or lines that can show strong chromospheric emission (such as He I D₃, the Na I D doublet, Hα, or the Ca infrared triplet) as well as regions with strong telluric features to prevent biases in the polynomial fit. The selected pixels were used to fit a second-order polynomial, and finally the template spectrum was divided by it.

3.2. Excluded regions

Since we use observations obtained from the ground to create the template, some regions of the spectra are affected by telluric contamination. Absorption and emission features coming from the Earth’s atmosphere on top of the stellar spectrum need to be excluded. The position of the telluric lines on the detector (pixel frame) is very stable. However, the position of the stellar spectral features on the detector depends on the relative velocity between the observer and the target, which is determined by the orbital and rotation motion of the Earth and by the absolute radial velocity of the target. The relative stellar velocity changes throughout the year, so the stellar spectrum shifts across the detector depending on the observation date. Except for targets with large radial velocity amplitudes (e.g. binaries), the main contribution to this velocity shift comes from the Earth orbital motion (barycentric Earth radial velocity, BERV), which can be up to approximately ±30 km s⁻¹, depending on the sky position of the target. To be sure any contribution from the Earth’s atmosphere at any moment of the year is removed, we disregarded wavelengths encompassing an enlarged region around each telluric feature and not only the telluric line itself. To define this region, we used the maximum BERV of the observations of the target used to build the template, which can be significantly smaller than 30 km s⁻¹, and we broadened each telluric line by subtracting and adding this maximum BERV to the blue and red line limits, respectively. In the VIS range, we masked the regions around telluric lines that are deeper than 5%. In the NIR range, the threshold was 4% in the Y and H bands and 2% in the J band, where the stellar lines are less abundant and shallower (see Nagel et al. 2019a, for details). Instead of completely masking telluric contaminated regions, we could attempt to correct shallow telluric features with techniques such as the one presented in Nagel et al. (2019a), which would increase the available spectral range. However, this is beyond the scope of this work and we leave it for a future paper. We also masked out regions close to chromospheric lines that can show a strong emission component in active stars, as was done before when normalizing the template orders.

3.3. Mask lines

3.3.1. Line identification

We used the normalised spectrum template to identify absorption features. Firstly, we searched for all local minima present in the template (except in the excluded regions mentioned above), which indicate the centre of the absorption lines. We considered the region around these minima (the pixels in the region going from the adjacent maxima at each side of a minimum, see the grey shaded region in Fig. 1) to represent the absorption lines of the stellar spectrum. Most of these spectral features do not come from a single elemental transition, but they are the result of line blending.

Fig. 1. Minima characterisation. The black dots and solid line are the template data showing an absorption feature. The shaded light grey region defines the pixels corresponding to the feature shown, which go from the local maxima (dotted grey vertical lines) at each side of the minimum (dotted dashed grey vertical line). The solid red line is the best Gaussian fit, and the dashed red vertical line indicates its minimum.

3.3.2. Line characterisation

For each identified minimum or line, we measured a set of parameters in the template to quantitatively characterise its shape. To have a more reliable estimation of these parameters, we fitted a Gaussian function to each spectral minimum

$$\begin{array}{c}\hfill G(\lambda )=d+aexp(-\frac{{(\lambda -b)}^2}{2c^2})·\end{array}$$(2)

Figure 1 shows an example of a line in a spectrum template with the best-fit Gaussian function. The pixels used in the fit are the ones between the two maxima at each side of the absorption feature.

To measure the central wavelength of each line, we used the position of the minimum of the Gaussian fit b. This gives a more robust and precise estimate than the wavelength corresponding to the local minimum of the spectrum template, which depends on the spectral sampling. The depth of each line is given by the flux at the Gaussian minimum, a + d, so small values correspond to deep lines. To measure the width of each line, we used the FWHM of the Gaussian fit, $2c\sqrt{2ln(2)}$. We computed the contrast or prominence of the line as the difference between the flux at the edges of the line (the adjacent local maxima) and the flux at the minimum of the line (what we defined as the line depth). This gives two values for the contrast; one corresponds to each side of the minimum, whose difference can be used to study the asymmetry of the line.

3.3.3. Line selection

Figure 2 shows the distribution of the parameters we defined above for the lines identified in the VIS template of the slowly-rotating (v sin i  ≤  2 km s⁻¹), M3.5 V, Luyten’s star (GJ 273, Karmn J07274+052). We obtained similar values for low v sin i targets across all spectral subtypes. As mentioned above, we excluded minima close to very strong lines and those that fall in regions with strong telluric contamination. We also excluded minima that do not contain enough pixels to be properly characterised (due to noise or too weak lines) and minima with profiles that highly deviate from a Gaussian (mainly blends for which the Gaussian fit gives unphysical values).

Fig. 2. Line parameters (top: FWHM, middle: depth, and bottom: contrast) as a function of the line central wavelength (left) and the number of lines as a function of the parameters, where the values of each order have been stacked (right) for Luyten’s star (Karmn J07274+035, M3.5 V, low v sin i). In all panels, colours correspond to different spectral orders. Lines in order-overlap regions, which are wider in the blue, are duplicated. The shaded grey areas in the left panels correspond to regions contaminated by telluric features, broadened by ±30 km s−1. The black lines in the right panels correspond to the distribution of the lines in all orders.

Based on these parameters, we defined some threshold values to select the lines that will become part of the final mask. We tried to select as many lines as possible so that the CCF has a high S/N, allowing us to derive its parameters with the best precision. At the same time, we want to avoid weak lines that contribute more noise than signal.

We tested several masks created using different parameter combinations (depth, width, contrast, asymmetry) and different selection values for each parameter. To analyse the performance of the masks, we used them to compute the CCFs of the observed spectra of some targets with low activity levels and low v sin i, and we checked which mask yielded the smallest RV scatter. This analysis shows that using different combinations of parameters or slightly different selection values did not significantly change the masks or the CCF parameters obtained. The final parameters and values we adopted are shown in Table 1.

We selected lines with FWHM within a certain range. The FWHM values had to be increased for targets with larger v sin i, up to few hundreds of km s⁻¹ for the most extreme cases, the ones with v sin i ∼ 30−40 km s⁻¹. In doing this, we made sure that most of the lines that were included in the CCF have a similar width while avoiding the narrowest lines, which tend to correspond to weak noisy lines that were not removed in previous steps, and the broadest lines, due to strong, unsuitable lines or blends. The removed lines do not provide significant RV information and they tended to increase the RV scatter. To avoid lines that are too asymmetric and to further remove weak lines, we found that the simplest method was to select lines for which the contrast at both sides is larger than a minimum value.

We also selected lines depending on their depth. However, for this parameter we did not use a single threshold value. The spectra of M dwarfs do not show a clear continuum due to the large number of features present; so across the spectral range we observe large differences in the line depths, with much deeper lines in regions containing strong molecular bands. This makes it difficult to just use a single cut value for the whole spectral range. Since we are using a normalised high S/N template, the depth values of the lines range from 1, for the shallowest lines, to 0 for the deepest ones. In some orders, all of the lines can be very shallow, with the deepest ones reaching values of only 0.8, for example, while other orders contain very deep lines, reaching values of 0.2, for instance. To minimise this bias, we rescaled the depth values so that, for each order, we considered that the depth goes from 1 to the average of the 10% deepest lines of the order. On this new scale, we selected lines with depths smaller than the 60% quantile of the new depth range.

We found that values around the ones shown in Table 1 allow us to select enough lines over the whole spectral range to obtain a CCF with a good S/N and precise parameters. Also, masks with these parameters worked well for most of the M dwarfs of the CARMENES survey sample. However, since the lines in the mask are weighted according to their RV information, changing these values to include more or less weak or wide lines did not change the CCF profiles significantly because these lines have lower weights (see next section).

Moreover, it is possible to fine-tune these values to a specific target and select optimum lines to study the velocity or the activity signals present in the spectrum. It has been shown that it is possible to classify spectral lines according to their sensitivity to stellar activity and, hence, according to their contribution to the activity signal in the final RVs (Wise et al. 2018; Dumusque 2018; Lisogorskyi et al. 2019; Cretignier et al. 2020). Using masks with lines selected to be more or less affected by activity should result in CCF profile parameters that are more or less sensitive to stellar activity variations.

3.4. Final mask

3.4.1. Position and weight

The final mask was comprised of the lines we selected following the approach explained above. Each line has a “position”, given by the minimum of the Gaussian fit in wavelength space, and a “weight”, which we calculated as the contrast divided by the FWHM of the line. In this way, the lines with the steepest profiles (i.e. the ones with more RV content) are the ones that contribute to the CCF the most.

When building the mask, all features are defined to have zero width (Dirac δ functions), independently of their measured values. At the time of computing the CCF of a specific observation, we forced each of the mask features to have the width of one pixel, which depend on the instrument and the position of the mask on the spectrum (see Sect. 4 for details).

In Fig. 3 we show a section of the normalised templates of stars of spectral types M0.0 V, M3.5 V, and M7.0 V. We also plot the local minima found and the final mask lines selected following the criteria explained here, along with their corresponding weight.

Fig. 3. Section of a normalised spectrum template used to build masks (solid black data points and line), minima found in the template (dashed blue line), and mask containing the minima selected based on some conditions (solid orange lines) for three stars of different spectral type and low v sin i: HD 238090 (M0.0 V; GJ 458A, Karmn J12123+544S; top panel), Luyten’s star (M3.5 V; GJ 273, Karmn J07274+052; middle panel), and Teegarden’s star (M7.0 V; Karmn J02530+168; bottom panel). The three targets are Doppler-shifted to the same reference frame.

3.4.2. Order merging

Since CARMENES is an echelle spectrograph, the observed spectra and also the high S/N template created from these observations are split in different spectral orders across the detector. The bluest orders have some overlap at the order ends, so there are some duplicated lines in our list. To have a single final mask covering the entire spectral range, we merged the lines in overlapping order regions, which occur in the VIS channel and in the bluest NIR orders. We identified lines that are the same in two consecutive overlapping orders by comparing their central wavelengths. If the lines have central wavelengths that are close enough to be the same line (we chose a maximum separation of 0.05 Å), we took the mean of their positions and their weights, and these are the values that are included in the mask. If a line in an overlapping region is not found in both orders that overlap, we did not include it in the final mask.

3.4.3. Absolute velocity

We used the SERVAL template made by adding observations to build the mask so the mask wavelengths obtained are shifted to the RV of the star. To remove this contribution, we could just Doppler-shift the line wavelengths by the opposite of the star’s velocity. But since these velocities are known to different accuracies, instead, we used a synthetic stellar spectrum as a reference. We cross-correlated the mask with a PHOENIX synthetic spectrum (Husser et al. 2013), which we assumed to define the zero-velocity frame, and we corrected the shift of the mask lines by the velocity obtained from the CCF. We used models with solar metallicity [Fe/H] = 0.0, surface gravity log g = 5.00, and effective temperature closest to the one of the star whose observations we used to build the template obtained from Schweitzer et al. (2019).

4. CCF and parameter computation

4.1. Spectrum preparation: Flux and wavelength correction

As before, we assumed that we have a high-resolution echelle spectrum, which is divided into different orders, and each order consists of a list of data points. One data point represents all of the detector pixels summed perpendicular to the dispersion direction, with a central wavelength and a flux value.

The CARMENES science spectrum is extracted relative to the spectrum of a flat lamp taken during the daily calibrations (Zechmeister et al. 2014). We divided the spectrum flux by an instrument response function (provided by the standard CARMENES reduction pipeline together with the reduced spectrum) to remove the SED of the flat lamp. Due to the way the spectrum was extracted, its flux values do not reflect the S/N of the original image. To weigh each order according to its original S/N, we rescaled the flux of each order so that its mean counts are equal to the square of the original S/N. When other spectral extraction pipelines are used (e.g. HARPS DRS, Mayor et al. 2003), the extracted spectral orders need to be corrected for the blaze response of the spectrograph grating. In addition, we applied a small correction to ensure a constant overall spectral energy distribution and thus minimise biases arising from different extinction values, for example (Berdiñas et al. 2016). These rescaling values were calculated from the spectrum with a higher S/N in the time series. By doing this, we altered the original S/N of the orders, but we avoided systematic biases in the CCFs.

In the bluest orders of the CARMENES VIS channel, the blue end of the order is noisier than the rest of the wavelengths. We found that removing the first 200−300 pixels improves the RV precision obtained.

We Doppler-shifted the spectrum to correct for the barycentric motion of the Earth and any instrumental drift measured so that the RV shift derived from the spectrum comes solely from the star. In doing so, the centre of the CCF is always located close to the absolute RV of the target.

4.2. Mask preparation: Line selection for spectra at different epochs

When computing the CCF of a series of observations of a target taken at different times of the year, we need to make sure that we are using exactly the same mask lines at each epoch to avoid measurement biases. Since the stellar spectrum shifts across the detector during the year according to the relative velocity between the observed star and the Earth, the lines at the detector edges are not always visible. Therefore, to make sure that we were using exactly the same mask lines at any epoch, we removed the mask lines with wavelengths smaller than the wavelength of the blue cutoff plus the maximum BERV of the target observations, and all the lines with a wavelength larger than the red cutoff minus the maximum BERV. Furthermore, since we displaced the mask across the spectrum to compute the CCF, some mask lines appear and disappear at the order cutoffs depending on the Doppler shift applied to the mask. So we also removed the lines that, depending on the mask position, do not always fall within the spectral order limits. In other words, for each order, we only kept the mask lines that are between λ_min and λ_max

$$\begin{array}{cc}\hfill {\lambda }_{\min }& ={\lambda }_1({\mathrm{BERV}}_{\max }-v_{\min })/c\hfill \\ \hfill {\lambda }_{\max }& ={\lambda }_n(-{\mathrm{BERV}}_{\max }-v_{\max })/c,\hfill \end{array}$$(3)

where λ₁ and λ_n are the minimum and maximum wavelength of the order, respectively, BERV_max is the maximum barycentric shift of the observations considered, v_min and v_max are the minimum and maximum shift applied to the mask (i.e. the CCF velocity range), and c is the speed of light.

When creating a mask from observations of a specific star, we already removed the stellar lines that fall on regions affected by tellurics at any time of the year. However, telluric lines overlap with different stellar lines depending on the absolute radial velocity of the star. One needs to further remove the mask lines that can be affected by tellurics if the absolute velocities of the star used to create the mask and the target star are different. We followed the same procedure as carried out when excluding regions during the mask creation process (Sect. 3.2, we broadened the telluric features by the maximum BERV of the target observations and removed the mask lines that overlap with them). Additionally, we needed to take into account the displacement of the mask across the spectrum. Some mask lines overlap with telluric-contaminated regions for some of the Doppler shifts applied, so we also removed them. After this process, we obtained a set of lines per order that, regardless of the Doppler shift applied, are not affected by tellurics and are visible in the spectrum at any epoch.

4.3. CCF location in RV space

Prior to computing the CCF of each order, we need to define the RV range and the steps that we use to Doppler-shift the mask, which determine the span and sampling of the CCF. We need to ensure the full coverage of the CCF span (its minimum and part of the wings), but also, we should not extend far into the continuum to save on computing time. To quickly locate the position of the CCF centre in RV space and estimate its width, we computed a test CCF spanning a large RV range of ∼200 km s⁻¹, using coarse steps of ∼1 km s⁻¹ in the spectral order with the highest S/N. When both the target’s absolute RV and v sin i are known, this first step can be skipped.

The sampling of the CCF velocity grid should be given by the average pixel size in velocity units, Δv, which can be obtained from the resolving power R of the instrument and the pixel sampling per spectral element. The resolving power at a given wavelength λ can be expressed as λ divided by the minimum difference that can be distinguished at this wavelength, δλ, which is typically the size of a spectral element. By applying the Doppler equation, we get the size of a spectral element in velocity units δv

$$\begin{array}{c}\hfill R=\frac{\lambda }{\delta \lambda }=\frac{c}{\delta v}\Rightarrow \delta v=\frac{c}{R}\Rightarrow \mathrm{\Delta }v=\frac{c}{R·s},\end{array}$$(4)

where c is the speed of light. By dividing this quantity by the average pixel sampling per spectral element, s, we get the size of a pixel in velocity units Δv, which gives us the step size by which we have to shift the mask. The average values for the CARMENES VIS and NIR arms are shown in Table 2. Another way to find this step size is to directly compute the average width of each spectrum pixel in velocity units. It is important to use the correct sampling to ensure that the CCF flux reflects the actual number of photons of the spectrum and, therefore, its S/N. This is used to estimate the statistical uncertainty on the RV. For example, if we were using a denser sampling (smaller RV steps), we would count the flux in the same pixels more than once, resulting in biased RV errors.

By using this RV step, we obtain a well-sampled CCF for which a Gaussian fit provides its minimum, width and depth, with their corresponding uncertainties. However, this type of sampling is not sufficiently dense to measure the bisector correctly and compute the BIS. For this, we used a denser sampling of ∼250 m s⁻¹. As with the RV, in order to estimate the uncertainty of the BIS, we needed to use the CCF computed with the correct sampling so that it is consistent with the original S/N of the spectrum (see Sect. 4.6.3 for details).

4.4. CCF computation details

The calculation of the CCF as a function of the Doppler-shift v was performed by applying Eq. (1) to each of the spectral orders separately. The product of the two signals, that is, the spectrum and mask, was computed in pixel space, so both the spectrum and mask have to be expressed in pixel units. Since we did not rebin or merge the spectral orders, each data point (wavelength and flux) corresponds to one pixel, so the spectrum is already in pixel space. The values of each data point correspond to the centre of the pixel. Our mask lines are Dirac δ functions (zero width) expressed in wavelength, so we mapped each one to the corresponding position in the spectral order, interpolating between the wavelengths of the spectral pixels. We then gave each of them a width of 1 pixel by adding or subtracting 0.5 pixels to the central value to obtain the pixel boundaries. In this way, we transformed our δ mask lines into bins of 1 pixel width.

As expressed in Eq. (1), if a spectrum pixel overlaps (totally or partially) with a mask line (non-zero region of the mask), we multiply the value of the flux in that pixel by the value of the mask line weight, thus taking the fraction of the pixel overlapping by the mask into account. If the overlap is not exact, we do not take the flux value of the centre of the pixel. Instead, we would linearly interpolate the flux between consecutive pixels, and take the interpolated flux value at the centre of the overlap region between the pixel and the mask. In this way, we would obtain a smoother CCF profile and still maintain the observed number of photons.

We computed the CCF on spectral orders that have been corrected for the blaze function of the echelle spectrograph to avoid biases due to the blaze slope. However, this does not take into account the difference in S/Ns along the order, from lower S/N edges to a peak S/N at the centre. Thus, a blaze-corrected spectrum does not preserve the flux and, therefore, its CCF does not immediately reflect the original S/N of the spectrum, which we need in order to estimate the uncertainty of the RV estimate. To correct for this, we should multiply all the pixels of an spectral absorption line by the value of the blaze function at the centre of the line. In this way, we approximately recover the original flux of the line without reintroducing the biasing effects caused by the slope of the blaze. However, since we did not know exactly which pixels correspond to a line in the spectrum, we instead multiplied the mask lines by the corresponding blaze value.

4.5. Final co-added CCF

We co-added the CCFs of the individual orders to obtain a final CCF that represents the full spectrum. Since we computed the order CCFs using the same RV grid, we could directly sum the flux at each RV step

$$\begin{array}{c}\hfill \mathrm{CCF}(v)={\displaystyle \sum _o{\mathrm{CCF}}_o(v),}\end{array}$$(5)

where CCF_o(v) is the CCF of order o and CCF(v) is the final co-added CCF.

Depending on the instrument and the star observed, not all the orders contain the same number of useable lines or have a sufficient S/N. For the CARMENES observations, some orders are severely affected by tellurics, which is reflected in a low number of lines in the mask and poor CCFs. Furthermore, for very late M dwarfs, the bluest orders show a markedly degraded S/N. The CCFs of very low S/N orders (S/N ≤ 10) or orders with few usable lines (≤10 lines) do not significantly contribute to the final CCF and may actually be counterproductive. Our tests show that discarding such poor orders may improve the profile of the final CCF, in the sense of obtaining better precision on the derived parameters.

We did not need to apply any weight to the CCFs of the different orders because the flux level of each of them is already proportional to the S/N of the original spectrum order, the number of lines used, and their quality, which are the parameters that we would use to asses the quality and, therefore, weight of each CCF. However, the final flux also reflects the values of the mask line weights. These weights depend on the RV-content of each line, but the absolute values are arbitrary and should not affect the final flux of the CCF. Therefore, we applied a normalisation. We only considered the mask lines that actually contributed to the co-added CCF, that is, the lines kept after removing telluric-contaminated regions and order ends of the orders used in the co-adding. We computed the average of the weights of these mask lines and divided the flux of the co-added CCF by it. In this way, the flux of the CCF reflects the true flux of the pixels used (scaled according to their RV-content) and the absolute weight of the mask lines is irrelevant.

4.6. CCF parameters and errors

To study the CCF profile, we fitted a Gaussian function to the central region (from maximum to maximum; see Fig. 4) of the co-added CCF. From this function, we derived the RV shift of the spectrum with respect to the mask, the FWHM, and the contrast of the CCF profile. We also analysed its asymmetry by studying its bisector. Here, we describe the estimation of these four parameters and their uncertainties.

Fig. 4. Illustration of the CCF activity parameters measured in this work. We plotted a typical CCF of an M dwarf (black dots) and the best Gaussian fit (red dashed line) together with the FWHM and contrast derived from the fit (left). We also plotted the bisector of the CCF (grey diamonds) and a zoom to the central region of the CCF to show the bisector shape (right). We indicate the top and bottom regions used to compute the BIS.

4.6.1. CCF minimum: RV

We estimated the velocity difference between the observed spectrum and the mask as the position of the minimum of the best Gaussian fit to the CCF. For time series RV work, we are interested in the relative velocities and, therefore, the differential values between the target and the mask are already appropriate. An intuitive method in order to estimate the RV uncertainty would be to use the formal error of the fit. However, the CCF is not a perfect Gaussian function as it represents the convolution of thousands of lines with different profiles. Therefore, the formal error of the fit overestimates the true uncertainty as the Gaussian model is not adequate. To obtain a realistic value for the uncertainty, we consider the total RV content in the CCF assuming photon noise. The capability of the CCF in providing a precise measure of the RV depends on its steepness and depth. A high S/N, deep, and narrow CCF is much more constraining than a low S/N, shallow, and wide CCF. In turn, this is defined by the S/N of the original spectrum by the number of lines used and by their profiles. Descriptions of the estimation of photon-limited errors in RV measurements from spectral lines or CCFs can be found in Butler et al. (1996), Bouchy et al. (2001), and Boisse et al. (2010), for example. Here we detail the procedure we followed, which is based on the methodologies of the previous references.

The first step is to compute the uncertainty of the flux on each sampling point of the CCF of each order. This uncertainty is propagated from the uncertainties of the individual pixels that are added when computing the CCF. We consider that the uncertainty on the flux of a pixel x in the spectral order o is given by its photon noise σ_PN, x and readout noise RON as follows:

$$\begin{array}{c}\hfill \sigma {{}_f}_{\mathrm{ox}}^2={({\sigma }_{\mathrm{PN},x})}^2+{\mathrm{RON}}^2=f_{\mathrm{ox}}+{\mathrm{RON}}^2.\end{array}$$(6)

Assuming Poissonian statistics, the photon noise of the pixel x is calculated as the square root of the total flux in the pixel, f_ox, where we use the flux value before rescaling the orders to a constant overall spectral energy distribution, that is, the value representing the observed S/N.

The uncertainty of the CCF of order o at velocity v, CCF_o(v), is then the square root of the quadratic sum of the uncertainties of all the pixels used when computing the CCF at that velocity. For partial pixels, the flux is weighted by the fraction of pixel that overlaps with the line l, which is similar to Eq. (1)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{CCF},o}^2(v)={\displaystyle \sum _{l=1}^m\sum _{x=1}^n{{\sigma }_f}_{\mathrm{ox}}^2{\mathrm{\Delta }}_{xl}(v).}\end{array}$$(7)

When computing the CCF of spectra extracted with a blaze function, the flux of pixel x is also normalised to the original flux level by multiplying it by the interpolated blaze value, corresponding to the pixel of the mask line l, so that the photon noise reflects the original S/N of the data. This is the same process as with the spectrum (see Sect. 4.4).

The uncertainty of each RV sampling point of the co-added CCF, ${\sigma }_{\mathrm{CCF}}^2(\mathit{v})$, is given by the quadratic sum of the uncertainties of all the orders used

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{CCF}}^2(v)={\displaystyle \sum _o{\sigma }_{\mathrm{CCF},o}^2(v),}\end{array}$$(8)

which is the error propagation of Eq. (5).

From the value of the uncertainty of the CCF flux (σ_CCF(v), the error on the y-axis), we estimate the velocity uncertainty of each CCF point (σ(v), the x-axis of the CCF) by using the derivative of the CCF as

$$\begin{array}{c}\hfill \sigma (v)={\sigma }_{\mathrm{CCF}}^2(v){\left(\frac{\mathrm{d}\mathrm{CCF}(v)}{\mathrm{d}v}\right)}^{-1}.\end{array}$$(9)

The uncertainty of the final RV is then given by the velocity errors of all the points of the CCF

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{RV}}={{\displaystyle \left(\sqrt{\sum _v\frac{1}{\sigma {\left(v\right)}^2}}\right)}}^{-1}·\end{array}$$(10)

For this procedure to provide realistic uncertainties, it is essential to employ the CCF sampled in RV steps corresponding to the average pixel size.

4.6.2. CCF profile variations: FWHM and contrast

To determine the width of the CCF, we used the FWHM of the best-fit Gaussian, which is measured in velocity units (typically km s⁻¹). The contrast of the CCF is the amplitude of the Gaussian divided by its baseline level. We inverted the sign of the amplitude, which is negative, and multiplied it by 100 so that the contrast values are positive and given in a percentage. We considered the formal errors of the fit as the uncertainties for these two quantities.

4.6.3. CCF profile variations: Bisector (asymmetry)

To quantify the asymmetry in the CCF, we used the standard bisector inverse slope as defined in Queloz et al. (2001). We measured the bisector of the CCF in the same interval in which we fitted the Gaussian, that is, the range comprised within the maxima at each side of the global minimum. We linearly interpolated the CCF profile and took the mid points between the left and right side profiles at about 100 different evenly spaced depths. This yields the bisector as a function of the depth d

$$\begin{array}{c}\hfill \mathrm{Bisector}(d)=(v({\mathrm{CCF}}_{\mathrm{left}}=d)+v({\mathrm{CCF}}_{\mathrm{right}}=d))/2,\end{array}$$(11)

where v(CCF_left = d) and v(CCF_right = d) are the interpolated velocities at depth d at the left and right sides of the CCF. We then consider the following two regions: one at the top of the CCF, from 90 to 60% of the baseline level, and one at the bottom, from 40 to 10% (see orange regions in Fig. 4). We computed the average RV of the points in these two regions, and the difference between them, RV_top and RV_bot, yields the value of the BIS

$$\begin{array}{cc}\hfill \mathrm{BIS}& ={\mathrm{RV}}_{\mathrm{top}}-{\mathrm{RV}}_{\mathrm{bot}}\hfill \\ \hfill & ={⟨\mathrm{Bisector}(d)⟩}_{d=60-90\%}-{⟨\mathrm{Bisector}(d)⟩}_{d=10-40\%}.\hfill \end{array}$$(12)

As mentioned above, using the CCF computed with the correct RV sampling (i.e. where the RV steps are the size of the average pixel width in velocity units) does not allow us to correctly measure the bisector. Instead, we used a CCF computed with a denser sampling, with steps of 250 m s⁻¹.

Regarding the uncertainty, since the BIS is computed from the (interpolated) points of the CCF, we started with the velocity uncertainties of each point of the CCF as computed above (Sect. 4.6.1). As was done with the CCF flux when computing the bisector, we linearly interpolated the RV errors of the actual CCF data points to the ∼100 points, corresponding to the evenly spaced depths used to define the bisector. For each bisector point, the error σ_bisector(d) comes from the interpolated uncertainties at each side of the CCF by error propagation of Eq. (11). We then computed the uncertainty on the top and bottom regions of the bisector as the mean of all the points considered in each region divided by the squared root of the number of points n_top or n_bot

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{top}}& ={⟨{\sigma }_{\mathrm{Bisector}}(d)⟩}_{d=60-90\%}/\sqrt{n_{\mathrm{top}}}\hfill \\ \hfill {\sigma }_{\mathrm{bot}}& ={⟨{\sigma }_{\mathrm{Bisector}}(d)⟩}_{d=10-40\%}/\sqrt{n_{\mathrm{bot}}}.\hfill \end{array}$$(13)

Finally the BIS uncertainty is the square root of the quadratic sum of the uncertainties on the top and bottom regions

$$\begin{array}{c}\hfill \sigma (\mathrm{BIS})=\sqrt{{\sigma }_{\mathrm{top}}^2+{\sigma }_{\mathrm{bot}}^2}.\end{array}$$(14)

5. Application to the CARMENES survey

5.1. Sample

We followed the procedure explained in the previous sections to build a set of masks using CARMENES survey observations to compute the CCFs for most of the targets in the sample. The CARMENES survey sample (Alonso-Floriano et al. 2015; Cortés-Contreras et al. 2017; Schweitzer et al. 2019; Reiners et al. 2018a) includes about 350 M dwarfs of subtypes from M0.0 to M9.0. Most of the targets are early and mid M dwarfs, and only less than 20 stars have a spectral type equal or later than M6.0. Most of the stars have low projected rotational velocities, v sin i ≤ 10 km s⁻¹, although some show v sin i values of a few tens of km s⁻¹. The sample includes both inactive and active stars for early and mid spectral types, while most of the late spectral type targets show high levels of activity.

In total, we computed CCF profiles and parameters of 323 targets in the VIS, excluding those with v sin i ≳ 60 km s⁻¹ and spectroscopic binary systems (already studied in Baroch et al. 2018). In the NIR, we computed CCFs for 305 targets. Since there is less RV content in the NIR than in the VIS, we could not get reliable masks and CCFs for the faintest targets and the ones with v sin i ≳ 25 km s⁻¹; we therefore excluded them from this analysis. We used different masks for targets of different spectral types, but this is especially the case for targets of different v sin i. In the following, we explain the tests we performed to select the final set of masks used for different types of targets, which we divide into slow and fast rotators.

5.2. Slowly-rotating stars

For each spectral subtype (M0.0 V to M0.5 V, M1.0 V to M1.5 V, etc.), we selected the objects with the highest S/N to build the best possible template. This generally corresponds to the brightest objects with the highest number of observations. We restricted this search to targets with low v sin i (v sin i ∼ 2 km s⁻¹) and low levels of activity. We did this because our early tests show that for targets with larger v sin i, we obtained better results (CCFs with smoother profiles and smaller RV dispersions) when using masks built from observations of a star with also larger rotational velocity, rather than using one of the standard masks obtained with observations of a slowly-rotating star. The masks were built using the selection parameters shown in Table 1.

Table 3 lists the targets used to create the masks together with the number of co-added observations, the S/N of the template, and the total number of lines of the masks for the CARMENES VIS and NIR spectral ranges. Figure 5 shows a histogram of the number of lines of the masks listed in Table 3 as a function of wavelength.

Fig. 5. Histogram of the line wavelengths of the masks used to compute the CCFs of the low v sin i targets for the VIS (top) and NIR (bottom) channels. Different colours correspond to the different spectral types of the targets used to build the masks (see targets and number of lines in Table 3). Light grey areas correspond to regions affected by telluric features, which are broadened by ±30 km s−1. Dark grey hatched areas correspond to regions without spectral information in the CARMENES data, corresponding to gaps between consecutive spectral orders (redwards of 11 500 Å), and intra-order gaps in the NIR detector mosaic.

We analysed the performance of the different VIS masks on targets of different spectral types. We did this for several targets with low v sin i, low activity levels, and more than 100 observations when possible, since most late spectral type targets are very active and have not been observed that many times. We computed the CCF of these targets with the different masks and studied the CCF profiles and the parameters time series.

Although they have a different number of lines, we find that for most of the early and mid spectral type stars, from M0 to M6–M7, masks built from early and mid spectral type targets result in smooth CCF profiles and RV time series with similar scatter and modulation. Therefore, we are able to obtain reliable CCF profiles and parameters with a single mask for all early and mid spectral type targets. However, we find that using more than one mask for different spectral subtypes slightly improves the precision of the parameters. We note that different masks result in RV time series with different absolute offsets, with differences of the order of a few tens of m s⁻¹, thus showing the difficulty in computing the absolute displacement of a mask.

Regarding the other CCF parameters, they also show similar dispersions and modulations for masks of different spectral types, and their absolute values also seem to depend on the mask used. FWHM mean values can vary by some hundreds of m s⁻¹ for masks of a different spectral type showing very similar RV scatters. The same happens with the contrast and the bisector. Contrast mean values can vary by about 10% and bisectors can fluctuate by some tens of m s⁻¹. This shows that the same mask must be used if we want to compare CCF parameters.

For early M dwarfs (M0 V−M1 V), the VIS channel CCF has the shape of an absorption line and rather flat wings. The same shape is found in G and K stars (see e.g. Fig. 1 in Mayor & Queloz 1995). As we move towards later-type stars, the flux level of the CCF continuum decreases and two humps appear at either side of the CCF (such as in Fig. 4). The appearance of the two humps and the decreasing continuum are a consequence of the very large density of spectral features in cooler stars. The mask lines are chosen to represent those spectral features that are best defined, that is, showing no obvious blends and asymmetries as well as with a clearer continuum. As the mask slides through the spectrum at different wavelength positions, the density of stellar lines is so high that a number of the mask lines randomly coincide with stellar absorption features. This means that the CCF continuum is depressed. The first CCF flux maximum occurs just before the mask lines enter the stellar absorption lines because all mask lines receive flux from the continuum as the mask is only composed of lines that are well behaved. The second maximum occurs for the same reason just after the mask lines leave the stellar lines. This explains the characteristic shape of the CCF for stars of spectral types later than M1 in the case of the CARMENES VIS spectra. Interestingly, we do not observe this phenomenon in the NIR channel CCFs as the lines in the spectrum are more widely spaced.

5.3. Fast-rotating stars

The spectra of stars with high v sin i values are quite different from those of low v sin i stars of similar types because their lines are much broader and most of them become heavily blended. Moreover, line profiles can be affected by distortions due to stellar activity, which is generally stronger in fast rotators. Using one of the “standard” masks made from a low v sin i and inactive target for a star with larger v sin i (from ∼3 km s⁻¹ up to ∼10 − 15 km s⁻¹) still delivers a smooth CCF profile. However, the RV dispersion is larger than if we were to use a mask made from an object of v sin i close to the one of the target star. This would be the case even if the spectral type is not a very good match.

For even faster rotators, using the standard masks results in a CCF profile that shows a global minimum close to the RV value of the target, but mostly this contains noise, with the CCF showing bumps across its profile. The profile becomes noisier as the v sin i of the star increases up to a point where we are not able to distinguish any minima for the stars with the most extreme v sin i values (v sin i ≳ 40 km s⁻¹). In using a mask that is also made from a large v sin i target, the CCF profile becomes less noisy, but we cannot obtain smooth profiles as is the case for slower rotators. When building such a mask, in order to be able to select any line, the FWHM cuts need to be larger than in the low v sin i case. The resulting masks have fewer lines, only a few hundred, which also impedes obtaining a smooth CCF because of the lower S/N that results from only considering a small number of lines. Another problem is that most of the fast-rotating stars are the ones with the latest spectral types and, therefore, they are fainter and deliver poorer S/N observations. Also, in general, fewer observations are available. All of this implies that the templates have a low S/N. In some extreme cases, we smoothed the templates applying a convolution with a Gaussian profile in order to identify local minima due to bona-fide absorption lines instead of noise features. For stars with v sin i ≳ 50 − 60 km s⁻¹, it becomes unfeasible to identify a large number of lines following our methodology. In all circumstances, and because of the numerous blends, CCF parameters of fast-rotating stars are much less reliable than for slow rotators.

5.4. Final masks

In light of these results, we used a grid of masks to compute the CCFs of the CARMENES survey sample. For targets with v sin i ≤ 4 km s⁻¹, we used the standard masks listed in Table 3. The procedure is to employ the mask corresponding to the spectral type closest to that of the target. For faster rotators, we used a different set of masks made from targets with larger v sin i values. We divided the targets according to their v sin i in bins of approximately 10 km s⁻¹, up to v sin i ∼ 60 km s⁻¹ in the VIS, and v sin i ∼ 30 km s⁻¹ in the NIR. For each bin, we selected targets of different spectral types for which we can obtain templates with a sufficiently high S/N and build the masks. We finally obtain one to three masks per v sin i bin covering different spectral types, depending on target availability.

5.5. Comparison between CCF and template-matching RVs

To analyse the performance of the CCF method regarding the precision of the resulting RVs, we compared them with the velocities computed using the standard CARMENES RV pipeline SERVAL (Zechmeister et al. 2018). As mentioned above, SERVAL computes the RV of a series of observations by fitting a high S/N template made from the observations themselves. It matches the template to the observations order-by-order, and the final RV is a weighted average of the RVs of the individual orders. Since it uses a complete spectrum template, it extracts more information than the CCF method, which only uses a subset of the absorption lines.

By default, SERVAL builds the template from the observations of the star that is being analysed, but it is also possible to use a different template built from observations of another star. Using another template gives better results in the case of faint stars or stars with a small number of observations. For the CCFs, we only used masks that were created from high S/N templates, except for a few fast rotators where the templates obtained have an S/N lower than the standard masks for slow rotators. Here, we use the default SERVAL results where a template has been created for each star with its own observations. All of the analysed stars have been observed at least five times.

Aside from using templates of a different S/N, another difference between both methods is the use of different orders to compute the final RV. The RV depends on the spectral region used since there is a chromatic dependence. Also, the RV dispersion appears to be larger in the blue region of the spectral range, where the activity effects are expected to be stronger (see e.g. Tal-Or et al. 2018). Moreover, ground observations, such as the CARMENES ones, show regions contaminated with tellurics. Some orders are more affected than others, and they may introduce additional scatter due to unidentified telluric features. In using different sets of orders and in specially removing blue orders in the case of active stars may reduce the RV dispersion. SERVAL uses the VIS orders 108 to 67, except in some faint and late spectral type targets, where about the 20 to 30 bluest orders are omitted. In the CCF VIS data, we used orders 108 to 68. However, between 5 and 10 orders were removed because the masks do not have enough lines to provide a reliable CCF. This is mainly due to the presence of tellurics. We also removed the bluest orders in the case of faint and late-type targets. Using these regions allowed us to obtain better results. In the case of the NIR arm, SERVAL mainly uses about 20 orders, which are not contaminated by tellurics, in the Y and H bands. We chose to use the same orders in the NIR CCFs to avoid telluric contamination.

In Fig. 6 we plotted the distribution of the RV time series dispersions obtained with the CCF method and with SERVAL’s template matching for each of the survey stars analysed in both VIS and NIR channels. We measured the dispersion as the weighted standard deviation of the RVs, where we used the squared inverse of the RV uncertainties as weights. We corrected the RVs for an instrumental nightly zero-point drift (Trifonov et al. 2018; Tal-Or et al. 2019). These nightly zero-points were calculated using the SERVAL RVs of stars with small variability observed in each night. We used the SERVAL measurements instead of the CCF ones because they are more precise in general. Figure 7 shows the correlation of the RV dispersions obtained with the two methods (the same quantities as in Fig. 6) for the analysed stars.

Fig. 6. Time series RV dispersions (weighted standard deviations) obtained with the CCF (black) and with SERVAL (orange) in the VIS (left) and NIR (right) channels. The dashed lines show the median of each distribution.

Fig. 7. Comparison of the time series RV dispersions (weighted standard deviation) obtained with the CCF and SERVAL methods, in the VIS (left) and NIR (right), which are colour-coded with spectral type. The grey dashed line indicates the one-to-one relation.

In the VIS channel, the RV dispersions of all the stars show a similar distribution for RVs calculated with SERVAL and CCFs. Although, in general, the scatter is slightly larger in the latter case. This is what we expected because the template matching technique used by SERVAL takes all of the spectral regions into account, rather than the regions around a set of some selected lines as the CCF does. Thus SERVAL is able to extract more RV information from the spectra.

For the NIR data, we obtain much larger dispersions than SERVAL in general. This may be caused by the low number of mask lines used (see Table 3), which in turn is due to a large number of regions containing tellurics in the NIR range. Not having enough lines in the mask results in CCFs with a low S/N and less precise parameters. We note that the NIR CCF results are not as optimised as in the VIS case, which is the one we used to perform most of our tests due to its better precision compared to the NIR. A different choice of parameters, both in the mask creation (different line selection values) and the CCF computation (different pixel and order masking) processes may lead to improved results. In the same vein, we expect the RV precision to be improved by using telluric-corrected spectra as presented in Nagel et al. (2019a).

There are some cases in which the dispersion of the CCF RVs is smaller than for SERVAL. This occurs for stars where the template has a poorer S/N, that is, in the case of faint stars or stars with a small number of observations. In Fig. 8, we plotted the ratio of the RV dispersions obtained with SERVAL and with the CCFs as a function of the number of observations co-added to build the SERVAL template for each star. We see that for the VIS data, the ratio is larger (which means lower dispersions for the CCF RVs compared to SERVAL RVs) in the case of stars with less than 10−20 observations. Late spectral type targets, which are faint, also show a large ratio. In these cases, instead of using a template made from the observations of the target itself, which is the default in SERVAL, we used a high S/N template built from observations of a similar star. We obtain RVs that are more precise and that show smaller scatter compared to the results of the default template.

Fig. 8. Ratio of the SERVAL and CCF RV dispersions (weighted standard deviations) as a function of the number of observations co-added by SERVAL to build the template, in the VIS (left) and NIR (right) channels, which are colour-coded with the stellar spectral type. The cyan crosses show the moving average of the data computed over windows of 20 data points.

5.6. Activity analysis

To show an example of the CCF activity indicators (FWHM, contrast, BIS), in Figs. 9 and 10, we plotted the RV time series, together with the three activity indicators, for an active star in the CARMENES survey sample, YZ CMi (GJ 285, Karmn J07446+035), corresponding to observations obtained with the VIS and NIR channels, respectively. YZ CMi is an active mid M dwarf (spectral type M4.5 V; Reid et al. 1995) with a rotational velocity of v sin i ≃ 4 km s⁻¹ (Reiners et al. 2018a) and a rotation period of P_rot ≃ 2.78 d (Díez Alonso et al. 2019).

Fig. 9. Time series of the four parameters derived from the CCF for the star YZ CMi in the VIS channel. We show RV, FWHM, contrast, and BIS (from top to bottom) as a function of barycentric Julian Date (left), phase of the rotation period (middle), and also plot the correlation of the activity indicators with the RV (right). The data points are coloured according to the phase.

Fig. 10. Same as Fig. 9, but for the NIR channel. In the phase-folded plots (central panels), we also plotted the VIS data as grey squares, which are shifted to the mean of the corresponding NIR parameter, to compare the amplitude of the modulation.

It shows a large-amplitude modulation in the RV of a few hundreds of m s⁻¹ peak-to-peak in both channels due to stellar activity. Both VIS and NIR data show similar behaviour. The BIS is the indicator that shows the clearest correlation with the RVs. Both parameters show a similar modulation in the phase-folded plot but in opposite phase, and the correlation plot shows a strong anti-correlation, with a Pearson’s correlation coefficient of about −0.9. This anti-correlation was first observed by Queloz et al. (2001) and has since been observed in several stars for which the RV modulation contains a signal caused by stellar activity (e.g. Bonfils et al. 2013, reported correlations in several M dwarfs). We also observed such anti-correlations in other active M dwarf stars from the CARMENES sample. The behaviour of the BIS is similar to the one observed for the chromatic index of YZ CMi, which is a measure of the variation of the RV with the wavelength (Zechmeister et al. 2018; Tal-Or et al. 2018).

The FWHM and the contrast also show a modulation at the stellar rotation period; although, it is not as clear as the BIS one. The correlation of these two parameters with RV is not very clear either, but there is some loop shape in which the data points follow the stellar rotation phase. This was already observed for another activity indicator, the differential line width, which is a parameter that measures variations in the stellar line width compared with a stellar template (again see Zechmeister et al. 2018).

This example shows the value of the CCF parameters to be reliable activity indicators. These activity proxies, together with other photospheric and chromospheric indicators of activity, are being used in several projects that include CARMENES observations to study stellar activity and to help distinguish between planetary and activity-related signals in RV time series. These projects include both M-dwarf data taken under the CARMENES survey program (Kaminski et al. 2018; Luque et al. 2018; Reiners et al. 2018b; Lalitha et al. 2019; Nagel et al. 2019b; Perger et al. 2019; Schöfer et al. 2019) as well as open time observations of other cool stars (Luque et al. 2019; Pallé et al. 2019). We observe that the indicators in different stars display notable differences depending on the spectral type and activity level of the stars. A study of such variation is, however, out of the scope of this paper, and it will be discussed in forthcoming articles.

5.7. Absolute RVs

We estimated absolute RVs of the CARMENES survey targets by using the RV values derived from the CCFs. Absolute RVs are necessary to determine stellar space motions and to carry out studies of galactic kinematics and dynamics. Since our template masks are calibrated against PHOENIX spectra, our CCF RVs are appropriate to determine the true offset of the wavelengths from their reference values. The radial velocity that we measured from the wavelength shifts of the absorption lines is displaced from the true motion of the stellar barycentre due to several effects (Lindegren & Dravins 2003). For cool dwarfs, the dominant effects are the gravitational redshift that the photons experience as they escape the potential well, which is created by the mass of the star, and the convective shift due to the motion of hot and cold material in the stellar photosphere. Therefore, apart from the RV that was directly derived from the CCF measurements, we also need to take these effects into account. In Appendix A, we describe the values used to compute the absolute RVs and the associated errors in detail, and we show the values obtained in Table A.1.

6. Summary

We present a procedure to build weighted binary masks, compute CCFs of high resolution echelle spectra with these masks, and derive RVs and the following standard CCF activity indicators: FWHM, contrast, and BIS. In order to build the masks, we used a stellar spectrum template built from observations and we selected the mask lines based on the profile of the local minima present in it.

We applied these methods to the CARMENES survey sample, which is a set of more than 300 nearby M dwarfs. By using CARMENES observations of the brightest and most frequently observed targets, we created several masks covering the different spectral subtypes and rotational velocities of the stars in the sample. We then used these masks to compute CCFs of the whole sample. Our tests show that the CCF parameters depend on the mask used, and they tend to be more precise if the mask was built from observations of a star with similar characteristics as the target. We found that the most important parameter to select a mask was the rotational velocity; results were more sensitive to this than to the spectral subtype of the mask. For fast rotators, the absorption lines in the spectrum become severely blended, so the number of mask lines that we were able to select decreases from a few thousand of the lowest v sin i masks to only a few hundred, making the CCF and its parameters much less reliable. Moreover, fast rotators tend to be active stars, so their line profiles are affected by deformations due to activity features on the stellar surface, which add noise to the results.

We compared the RVs derived from the CCF method with the ones obtained with the standard CARMENES pipeline SERVAL, which performs a fit of the observation with a spectrum template obtained by co-adding the observations themselves. As expected, SERVAL generally performs better than the CCF because it is able to exploit more information from the whole spectrum, while the CCF only uses the region around a set of a few selected lines. However, in some cases, the scatter was smaller in the CCF RVs than in SERVAL. This happened in faint stars or stars with few observations, where the template employed by SERVAL has a low S/N, while the CCF uses masks built from high S/N templates.

We show the behaviour of the CCF activity indicators of YZ CMi, wich is an active mid-type M dwarf that displays a large modulation due to activity in its RV time series. We found that the BIS shows a clear anti-correlation with the RVs in a similar way as the chromatic index. On the other hand, the FWHM and contrast show a more complicated correlation with the RVs, but if we consider the stellar rotation phase of the signals, both parameters show phase-shifted correlations.

Finally, we estimated the absolute RVs of the CARMENES sample stars using the CCF RVs and discuss the main sources of shifts that affect the M-dwarf spectra. This includes the redshift caused by the gravitational potential of the star and the shift due to convective motions in the photosphere.

In conclusion, even though other algorithms allow us to obtain more precise RVs in the case of cool stars, the CCF method still offers a fast and straightforward way of obtaining RVs by simply cross-correlating the observed spectrum with a mask, which is a simple model of the absorption lines in the spectrum. Moreover, we can compute different indicators related to stellar activity from the CCF profile, which have proven to be valuable when studying activity-related signals in RV time series.

Acknowledgments

We thank the anonymous referee for the constructive report. 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 FICTS-2011-02 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, Institut 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 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. Based on data from the CARMENES data archive at CAB (INTA-CSIC).

References

Appendix A: Absolute RVs

To compute the absolute RVs of the CARMENES survey stars, RV_abs, we followed

$$\begin{array}{c}\hfill {\mathrm{RV}}_{\text{abs}}={\mathrm{RV}}_{\mathrm{CCF}}+v_{\mathrm{grav}}+v_{\mathrm{conv}},\end{array}$$(A.1)

where RV_CCF is the RV directly derived from the CCF measurements, v_grav is the gravitational shift of the star, and v_conv is the convective shift of the stellar photosphere.

We computed RV_CCF as the weighted mean of the RV time series of each target. Ideally, we would add in the model of each time series signal from companions and stellar activity on top of the RV zero point. For simplicity, we ignored those effects, which are mostly short-term (and therefore averaged out in long time series) and at the m s⁻¹ level. Exceptions include single-lined spectroscopic binaries with long periods where the system velocity has not been established yet. Double-lined spectroscopic binaries in the sample were not considered here; Baroch et al. (2018) already provided absolute RVs for them. We used the nightly zero-point corrected VIS data and removed the observations without drift correction (except in stars without any instrumental drift estimation). We estimated the uncertainty arising from the dispersion of the different RV measurements of each star, σ_CCF, as their weighted standard deviation divided by the square root of the number of measurements.

A systematic error source when estimating the absolute RV shift is the zero-point of the CCF masks. As explained in Sect. 3.4.3, in order to transform the masks to an absolute reference, we cross-correlated them with PHOENIX synthetic models, which we assumed to have a zero offset, and we also used the RV derived from the best-fit Gaussian to the CCF to correct the mask shift. There is uncertainty as to the mask zero-point, σ_mask, because of mismatches between the mask and the model lines. To measure this uncertainty, we used the formal error of the Gaussian fit, which has typical values of ∼30 m s⁻¹ for a mask made from low v sin i targets, and it increases by some tens of m s⁻¹ for the fastest rotators.

The gravitational redshift of a photon escaping from the surface of a star of mass M and radius R is given by

$$\begin{array}{c}\hfill v_{\mathrm{grav}}=\frac{\mathit{GM}}{\mathit{cR}},\end{array}$$(A.2)

where G is the gravitational constant and c is the speed of light. We used this equation with stellar mass and radius values from Schweitzer et al. (2019) to estimate the gravitational redshift, and we assumed the uncertainties derived by error propagation from Eq. (A.2).

The net velocity shift due to convective motions in the photosphere is more difficult to quantify. In the Sun, the effect of the upward motion of the granules and the downward motion of the intergranular lanes results in a net average blueshift of about 400 m s⁻¹. However, the exact value depends on the characteristics of the line used, such as the line strength, the excitation potential, or the formation depth in the atmosphere. There is also a limb effect: The value of the net shift changes from the centre to the limb of the solar disc (Löhner-Böttcher et al. 2019). For other cool stars, 3D hydrodynamical atmosphere models that compute the convective shift using several lines in a specific wavelength range (Allende Prieto et al. 2013), as well as observations of specific spectral lines (Meunier et al. 2017), also show a dependence of the convective shift on the line or lines used. In particular, lines at shorter wavelengths experience larger blueshifts than those at longer wavelengths, and shallower lines (formed deeper in the atmosphere) are also more blueshifted than stronger lines. The effective temperature of the star also affects the convective shift (Bauer et al. 2018) from blueshift values of 500 m s⁻¹ for early G-type stars decreasing until values of 100 m s⁻¹ for late K-type stars. The blueshift also decreases with the metallicity, the surface gravity, and the activity level of the star.

We found neither theoretical nor observational studies that delved into the effects of convective blueshift in M dwarfs. Therefore, we extrapolated the results obtained for earlier spectral type stars. We expect the convective shift in M dwarfs to be small since they have lower effective temperatures and can have large levels of activity. There even seems to be evidence of a reversed convection pattern, showing a convective redshift instead of a blueshift (Kürster et al. 2003). Due to the difficulty in estimating a value for the convective shift v_conv for our target stars, we assumed it to be zero with an uncertainty σ_conv of 100 m s⁻¹.

The final absolute RV values, RV_abs, obtained using Eq. (A.1), are shown in Table A.1. The RV_abs uncertainties are the quadratic sum of the following four sources of error that we considered: the error of the RV measurements dispersion σ_CCF, the error from the gravitational redshift σ_grav, 100 m s⁻¹ to account for the uncorrected convective shift σ_conv (which for most targets is the dominating one), and the error of the zero-point of the mask σ_mask

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{abs}}^2={\sigma }_{\mathrm{CCF}}^2+{\sigma }_{\mathrm{grav}}^2+{\sigma }_{\mathrm{conv}}^2+{\sigma }_{\mathrm{mask}}^2.\end{array}$$(A.3)

The same table also contains the average RV from the CCF measurements, RV_CCF (which already includes the correction of the mask zero-point), together with the combined uncertainty of the RV scatter σ_CCF and the mask zero-point σ_mask given by their quadratic sum. We also list these two uncertainties separately. Finally, we also give the gravitational redshift estimates together with their uncertainties.

All Tables

All Figures

Fig. 1. Minima characterisation. The black dots and solid line are the template data showing an absorption feature. The shaded light grey region defines the pixels corresponding to the feature shown, which go from the local maxima (dotted grey vertical lines) at each side of the minimum (dotted dashed grey vertical line). The solid red line is the best Gaussian fit, and the dashed red vertical line indicates its minimum.

In the text

Fig. 2. Line parameters (top: FWHM, middle: depth, and bottom: contrast) as a function of the line central wavelength (left) and the number of lines as a function of the parameters, where the values of each order have been stacked (right) for Luyten’s star (Karmn J07274+035, M3.5 V, low v sin i). In all panels, colours correspond to different spectral orders. Lines in order-overlap regions, which are wider in the blue, are duplicated. The shaded grey areas in the left panels correspond to regions contaminated by telluric features, broadened by ±30 km s−1. The black lines in the right panels correspond to the distribution of the lines in all orders.

In the text

Fig. 3. Section of a normalised spectrum template used to build masks (solid black data points and line), minima found in the template (dashed blue line), and mask containing the minima selected based on some conditions (solid orange lines) for three stars of different spectral type and low v sin i: HD 238090 (M0.0 V; GJ 458A, Karmn J12123+544S; top panel), Luyten’s star (M3.5 V; GJ 273, Karmn J07274+052; middle panel), and Teegarden’s star (M7.0 V; Karmn J02530+168; bottom panel). The three targets are Doppler-shifted to the same reference frame.

In the text

Fig. 4. Illustration of the CCF activity parameters measured in this work. We plotted a typical CCF of an M dwarf (black dots) and the best Gaussian fit (red dashed line) together with the FWHM and contrast derived from the fit (left). We also plotted the bisector of the CCF (grey diamonds) and a zoom to the central region of the CCF to show the bisector shape (right). We indicate the top and bottom regions used to compute the BIS.

In the text

Fig. 5. Histogram of the line wavelengths of the masks used to compute the CCFs of the low v sin i targets for the VIS (top) and NIR (bottom) channels. Different colours correspond to the different spectral types of the targets used to build the masks (see targets and number of lines in Table 3). Light grey areas correspond to regions affected by telluric features, which are broadened by ±30 km s−1. Dark grey hatched areas correspond to regions without spectral information in the CARMENES data, corresponding to gaps between consecutive spectral orders (redwards of 11 500 Å), and intra-order gaps in the NIR detector mosaic.

In the text

	Fig. 6. Time series RV dispersions (weighted standard deviations) obtained with the CCF (black) and with SERVAL (orange) in the VIS (left) and NIR (right) channels. The dashed lines show the median of each distribution.
In the text

	Fig. 7. Comparison of the time series RV dispersions (weighted standard deviation) obtained with the CCF and SERVAL methods, in the VIS (left) and NIR (right), which are colour-coded with spectral type. The grey dashed line indicates the one-to-one relation.
In the text

	Fig. 8. Ratio of the SERVAL and CCF RV dispersions (weighted standard deviations) as a function of the number of observations co-added by SERVAL to build the template, in the VIS (left) and NIR (right) channels, which are colour-coded with the stellar spectral type. The cyan crosses show the moving average of the data computed over windows of 20 data points.
In the text

	Fig. 9. Time series of the four parameters derived from the CCF for the star YZ CMi in the VIS channel. We show RV, FWHM, contrast, and BIS (from top to bottom) as a function of barycentric Julian Date (left), phase of the rotation period (middle), and also plot the correlation of the activity indicators with the RV (right). The data points are coloured according to the phase.
In the text

	Fig. 10. Same as Fig. 9, but for the NIR channel. In the phase-folded plots (central panels), we also plotted the VIS data as grey squares, which are shifted to the mean of the corresponding NIR parameter, to compare the amplitude of the modulation.
In the text