Issue 
A&A
Volume 616, August 2018



Article Number  A155  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201731010  
Published online  07 September 2018 
The GAPS Programme with HARPSN at TNG
XVII. Line profile indicators and kernel regression as diagnostics of radialvelocity variations due to stellar activity in solarlike stars^{⋆,}^{⋆⋆}
^{1}
INAF – Osservatorio Astrofisico di Catania, Via S. Sofia 78„ 95123 Catania, Italy
email antonino.lanza@inaf.it
^{2}
Dipartimento di Fisica e Astronomia Galileo Galilei – Università di Padova, Via Francesco Marzolo 8, 35122, Padova, Italy
^{3}
INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122, Padova, Italy
^{4}
INAF – Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34143 Trieste, Italy
^{5}
INAF – Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, Italy
^{6}
INAF – Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, 80131 Napoli, Italy
^{7}
European Southern Observatory, Alonso de Córdova 3107, Santiago, Vitacura, Chile
^{8}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150762 Porto, Portugal
^{9}
INAF – Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy
^{10}
Fundación Galileo Galilei  INAF, Rambla José Ana Fernandez Pérez 7, 38712 Breña Baja, TF, Spain
^{11}
INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy
^{12}
INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate (LC), Italy
^{13}
Instituto de Astrofísica de Canarias, C/ Via Láctea s/n, 38205 La Laguna, TF, Spain
^{14}
Departamento de Astrofísica, Univ. de La Laguna, Av. del Astrofísico Francisco Sánchez s/n, 38205 La Laguna TF, Spain
^{15}
Dipartimento di Fisica, Università di Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Roma, Italy
^{16}
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
^{17}
INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
Received:
20
April
2017
Accepted:
10
April
2018
Aims. Stellar activity is the ultimate source of radialvelocity (hereinafter RV) noise in the search for Earthmass planets orbiting latetype mainsequence stars. We analyse the performance of four different indicators and the chromospheric index log R^{′}_{HK} in detecting RV variations induced by stellar activity in 15 slowly rotating (υ sin i ≤ 5 km s^{−1}), weakly active (log R^{′}_{HK} ≤ −4.95) solarlike stars observed with the highresolution spectrograph High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPSN).
Methods. We consider indicators of the asymmetry of the crosscorrelation function (CCF) between the stellar spectrum and the binary weighted line mask used to compute the RV, that is the bisector inverse span (BIS), ΔV, and a new indicator V_{asy(mod)} together with the full width at half maximum (FWHM) of the CCF. We present methods to evaluate the uncertainties of the CCF indicators and apply a kernel regression (KR) between the RV, the time, and each of the indicators to study their capability of reproducing the RV variations induced by stellar activity.
Results. The considered indicators together with the KR prove to be useful to detect activityinduced RV variations in ~47 ± 18 percent of the stars over a twoyear time span when a significance (twosided pvalue) threshold of one percent is adopted. In those cases, KR reduces the standard deviation of the RV time series by a factor of approximately two. The BIS, the FWHM, and the newly introduced V_{asy(mod)} are the best indicators, being useful in 27 ± 13, 13 ± 9, and 13 ± 9 percent of the cases, respectively. The relatively limited performances of the activity indicators are related to the very low activity level and υ sin i of the considered stars. For the application of our approach to sunlike stars, a spectral resolution allowing λ/Δλ ≥ 10^{5} and highly stabilized spectrographs are recommended.
Key words: planetary systems / stars: activity / stars: latetype / starspots / stars: atmospheres / techniques: radial velocities
Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.
The IDL macro to compute the line profile indicators is also available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/616/A155
© ESO 2018
1. Introduction
The search for rocky planets around latetype stars has pushed the measurement of stellar radial velocity (hereafter RV) into the m s^{−1} regime because of the very low mass ratio of those planets relative to their host stars (e.g. Mayor et al. 2009; Queloz et al. 2009; Pepe et al. 2011; Díaz et al. 2016). At that level of precision, even the most stable stars show intrinsic RV variations on a variety of timescales, ranging from a few minutes, characteristic of pmode oscillations, to decades, associated with activity cycles produced by the modulation of their surface magnetic fields.
Variations on timescales from minutes to hours, induced by oscillations and photospheric convection (granulation), can be significantly reduced by averaging RV measurements collected with a suitable cadence (Dumusque et al. 2011a). The effects of photospheric active regions (hereafter ARs), with a lifespan comparable with the stellar rotation period or longer, are much more subtle and difficult to identify and remove. The magnetic fields of ARs produce a perturbation of the local brightness that induces distortions on the profiles of photospheric lines (the socalled flux effect) and an attenuation of the local convective blueshifts of the spectral lines, which manifests itself as an apparent redshift in the discintegrated RV (see Lagrange et al. 2010; Meunier et al. 2010; Lanza et al. 2011, and references therein). The flux effect is roughly proportional to the filling factor of the ARs and the projected rotation velocity of the star, that is its υ sin i (Desort et al. 2007), while its sign depends on the position on the stellar disc and the contrast of the brightness inhomogeneities, being opposite for cool spots and bright faculae. On the other hand, the convective shift effect has always the same sign for both kind of inhomogeneities.
In the Sun, photospheric faculae have a total area approximately one order of magnitude larger than sunspots, while their contrast is very low at the disc centre and increases towards the limb (cf. Unruh et al. 1999). Together with the small υ sin i, this makes the convective shift effect the dominant RV perturbation in the models computed by Lagrange et al. (2010) and Meunier et al. (2010). Their predictions have been confirmed by observations on timescales ranging from a few rotation periods to the elevenyear activity cycle. The amplitude of the RV variation of the Sun as a star can reach up to ~10−12 m s^{−1} on timescales significantly shorter than the elevenyear cycle, in some cases as short as a few tens or hundreds of days (Haywood et al. 2016; Lanza et al. 2016). On a timescale of one week, an rms of 1.33 m s^{−1} has been measured with the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPSN; Dumusque et al. 2015). Those amplitudes are comparable to or larger than those expected from an Earthmass planet orbiting an F, G, or Ktype star with a period of days, thus making stellar activity the most challenging source of false positives in the detection of telluric planets around solarlike stars (Fischer et al. 2016).
Several approaches have been proposed to mitigate the effects of activity on RV time series. They are particularly useful in the case of latetype stars whose ARs are stable for at least a few consecutive rotations, thus allowing us to constrain their longitudes and parameters from spectroscopic or simultaneous photometric data (e.g. Boisse et al. 2009; Lanza et al. 2011; Haywood et al. 2014; Donati et al. 2014), while they are of limited use for stars whose AR lifespan is shorter than the rotation period as happens in the Sun.
In this work, we shall focus on the use of indicators of line profile asymmetries as can be derived from the crosscorrelation function (CCF) between the stellar spectrum and a binary weighted line mask, specifically that used to measure the RV itself (cf. Sect. 2 for details on the computation of the CCF). All the RV information content of the spectrum is contained in the CCF. Fiberfed stabilized spectrographs, such as HARPS or the Spectrographe pour l’Observation des Phénomènes des Intérieurs stellaires et des Exoplanètes (SOPHIE), have a very stable wavelength reference frame and provide a CCF with a signaltonoise ratio typically ≥10^{3} (Udry et al. 2006). The subtle line profile asymmetries induced by stellar ARs in the case of a moderately rotating star (υ sin i ≲ 10−15 km s^{−1}) cannot be identified in the individual lines, while they can be detected in the high signaltonoise CCF.
The first indicator of line profile asymmetry was the line bisector (Queloz et al. 2001), which is a sensitive indicator in the case of rapidly rotating stars whose variability is dominated by cool spots. The amplitude of its variation is approximately proportional to the spot filling factor and the square of projected rotational velocity υ sin i (cf. Hébrard et al. 2014, and references therein). Therefore, for slowly rotating stars (υ sin i ≲ 5−7 km s^{−1}), other indicators have been introduced in an attempt to improve sensitivity to activityinduced variations (e.g. Boisse et al. 2011; Figueira et al. 2013). The main advantage over chromospheric indicators is that asymmetry indicators measure the variation of the same CCF profile used to derive the stellar RV. A CCF can be regarded as an average photospheric line which, thanks to its very high signaltonoise ratio, may provide a better correlation with the RV than the chromospheric proxies whose variations tend to level off at very low levels of activity. Moreover, the distribution of the ARs on the photosphere may sometimes differ from that in the chromosphere. In addition to CCF asymmetry indicators, the full width at half maximum (FWHM) of the CCF itself has provided a useful proxy for the activityinduced RV variations in slowly rotating stars (e.g. Dumusque 2014; Santos et al. 2014; Fischer et al. 2016), so we shall include it in our investigation. We stress that a measurement of the mean line profile distortions by means of the CCF variations would be the ideal activity indicator. It is the CCF variations that create the measured shift in RV, therefore by characterizing the former we can characterize the latter. This is completely different than associating the RV signal with chromospheric indicators or photometry, which are proxies of activity (i.e. created by the same physical cause that creates RV variations, that is surface magnetic fields) but which, being proxies, might not always be present or correlated with the RV variations.
We analyse a dataset of spectra of Gtype mainsequence stars with known exoplanets to further explore the advantages and the drawbacks of CCF line profile indicators. Our sample consists of stars with a low level of activity and slow rotation (2 ≲ υ sin i ≲ 5 km s^{−1}), comparable to the Sun, that have been observed with HARPSN at the Italian 3.58m Telescopio Nazionale Galileo within the coordinated observational programme Global Architecture of Planetary Systems (GAPS; Covino et al. 2013; Desidera et al. 2013). The aim of this specific project is the detection of additional planets in those systems. While our new discoveries have been reported in other works (e.g. Desidera et al. 2014; Damasso et al. 2015b), here we focus on the study of the correlations of different line profile asymmetry indicators and the FWHM of the CCF with the residual RV variations, that is, those obtained by subtracting the modulations due to the known planets, and whose origin can be confidently attributed to stellar activity.
2. Observations
Our sample consists of 15 latetype stars hosting planets listed in Table 1 together with their effective temperature T_{eff}, surface gravity log g, metallicity [Fe/H], projected rotational velocity υ sin i, mean chromospheric Ca II H&K index with its standard deviation, number N_{RV} of RV measurements considered for the present investigations (see below), and their total time span Δt (see below). The values of T_{eff}, log g, and [Fe/H] were extracted from the reference in the last column, while the values of the υ sin i came from the exoplanets.org database (Han et al. 2014), except for XO2S for which we used the parameters obtained by Damasso et al. (2015a) and Biazzo et al. (2015), and HD 108874 for which we used Benatti et al. (2017). Typical uncertainties in log g range between 0.06 and 0.12, while for [Fe/H] they range between 0.03 and 0.07. We used the spectra gathered within the GAPS programme to compute the mean value of and its standard deviation following the method described by Lovis et al. (2011).
Sample of stars for which we consider asymmetry indicators and FWHM of the CCF.
HARPSN is a fiberfed crossdisperser echelle spectrograph with a fixed setup covering the spectral range 383–690 nm with a resolution of 115 000 (Cosentino et al. 2012, 2014). We used exposure times of 600 or 900 s yielding a signaltonoise ratio (hereafter S/N) between ~80 and ~210 at 550 nm, according to the apparent brightness of the target. The data reduction software (hereafter DRS) of HARPSN is the pipeline used to obtain the RV by crosscorrelating the observed spectrum with a binary weighted line mask. These masks are chosen from templates matching the spectral types of the observed stars and the method used to compute the CCF is described in Baranne et al. (1996). The DRS provides also the CCF profile, which is used to compute the asymmetry indicators and the FWHM as explained in Sect. 3. The median error of the RVs as returned by the DRS ranges from 0.45 to 1.06 m s^{−1} for the stars in our sample, except for X02S for which it is 2.34 m s^{−1}.
HARPSN is enclosed in a controlled environment that ensures a stability of the RV measurements better than 1 m s^{−1} along several years. The uniform illumination of the entrance pupil and the stability of the system provide also a very stable CCF that can be used for our purposes. A very small variation of the focus of the instrument was detected during the first period of its operation. This variation perturbed the CCF by slightly modifying its FWHM without appreciably affecting the RV measurements. The focus was stabilized and the FWHM trend eliminated starting from JD 2456738. Therefore, we consider only spectra taken after that date with a time baseline of about two years. We discard a few measurements obtained at large airmass (≳1.8) and/or at small signaltonoise ratio (S/N ≲ 25−30 at 550 nm).
We subtracted the orbital motion of the known planets from the RV time series of each star, by performing an independent fit of the joined literature and HARPSN datasets. The number of planets in each system was taken from the literature with the addition of an independent RV offset and a jitter for each dataset, and of a linear longterm trend. We computed the RV residuals by choosing the median of the posterior distributions of the orbital parameters, computed with the PyORBIT^{1} code as described in Malavolta et al. (2016). The analysis of the full RV time series and the corresponding implications in terms of orbital parameters of known planets and detectability of additional companions will be presented in a forthcoming study. For the stars in our sample, the RV modulations induced by the known planets are at least an order of magnitude greater than the RV fluctuations produced by activity, thus the impact of the latter on the orbital parameters of the known planets is small. Nevertheless, a bias is introduced in the RV residuals when we fit a keplerian model because the Fourier components of the stellar activity variation at the orbital periods of the planets and their harmonics can affect the keplerian fit. To quantify this effect, we repeated the analysis by including in the model a linear correlation between the RV variation and the index of chromospheric activity , at least for the subsets for which this activity indicator is available, and compared the results with the keplerian model without this correlation (see Sect. 4).
The residuals of the model including the linear correlation between the RV and the chromospheric index were searched for periodicities longer than 1.0 d arising from possible additional planets by means of the Generalized LombScargle periodogram of Zechmeister et al. (2009). No significant periodicities with a falsealarm probability (FAP) < 0.001 were found. The FAP estimated with the analytical formula of Zechmeister et al. (2009) was checked by computing the periodograms for 10 000 random shufflings of the RV residuals, while the times of the observations were held fixed. Therefore, we assume that the residual RV time series are dominated by the intrinsic RV variation of the stars and use them to study their correlations with the activity indicators as explained in Sect. 3.
The standard deviations of the RV residuals of the models including only the keplerian motions and a longterm linear trend in the RV range from ~1.47 m s^{−1} for HD 99109 to ~5.73 m s^{−1} for HD 75898, the latter being one of our fastest rotators with υ sin i = 4.5 ± 0.5 km s^{−1} (cf. Table 1). Considering the models including also the correlation between the RV and the activity index, we see that this correction practically affects only the residuals of HD 75898 because their standard deviation is decreased to ~5.58 m s^{−1}. In the case of the Sun, the standard deviation of the RV variations ranges between 2.46 and 3.44 m s^{−1} on comparable timescales, depending on the level of activity (cf. Sect. 4.2 of Lanza et al. 2016), which is similar to the average of the present sample.
3. Methods
Several indicators of spectral line asymmetry have been proposed in the literature (e.g. Boisse et al. 2011; Figueira et al. 2013). The bisector inverse span (hereafter BIS) is probably the most widely used and proved to be well correlated with the activityinduced RV variations of active stars with υ sin i ≳ 7−10 km s^{−1} (e.g. Queloz et al. 2001). In stars with a smaller rotation velocity and a lower level of activity, the BIS is not always the best asymmetry indicator and other indicators have been explored (Figueira et al. 2013). Following those results, in addition to the BIS, we consider the two CCF asymmetry indicators ΔV of Nardetto et al. (2006) and V_{asy} of Figueira et al. (2013), the latter designed to fully exploit the information content of the different portions of the CCF (see Sect. 3.1 for their definition).
All the stars in our sample have υ sin i ≲ 5 km s^{−1} and a typical chromospheric index , that is, an activity level comparable with that of the Sun close to the minimum of the elevenyear cycle^{2}. The small filling factor of their active regions (≲1 percent) and the small rotational broadening of their line profiles in comparison to the width of the spectrograph instrumental profile make line distortions barely detectable, thus making all asymmetry indicators remarkably less sensitive than in more rapidly rotating and more active stars. In that case, the FWHM may become a more sensitive activity proxy, so we consider it in addition to the line asymmetry indicators. The better performance of the FWHM is likely related to the dominance of facular areas in lowactivity stars (e.g. Dumusque 2014; Dumusque et al. 2014). Faculae induce a variation of the convective shifts of the line profiles leading to a perturbation that is not localized in the core of the lines, but affects to some extent also their wings and therefore the FWHM of the CCF.
3.1. Asymmetry indicators and FWHM of the CCF
To compute the BIS, we use the CCF profile provided by the DRS and expressed as the crosscorrelation value per RV bin of width 250 m s^{−1}. We fit a Gaussian profile with a variable continuum level to the CCF to determine the continuum level itself and use it to normalize the CCF profile. The minimum of the normalized CCF profile gives the depth D of the CCF, also referred to as the CCF contrast. It is divided into 100 intervals and, for the flux level corresponding to each interval, the RV abscissa is quadratically interpolated on the red and the blue wings of the profile, respectively. The bisector at a given flux level is defined as the arithmetic mean of the RV corresponding to the flux level on the red and blue wings, respectively. Finally, the BIS is computed as the difference of the mean RV of the bisectors in the top and the lower parts of the CCF, where the top part consists of the elements between 10 and 40 percent of D below the continuum, while the lower part is taken between 60 and 90 percent of D below the continuum.
In principle, the error of the BIS can be computed by considering the maximum of two different estimates, that is the mean standard error of the 30 differences between the bisector in the top and the lower parts of the CCF and times the formal RV error as given by the DRS. This approach assumes that: a) all the 30 differences between the top and lower parts of the bisector are independent, which is not the case because the points along the CCF are correlated owing to the crosscorrelation procedure giving the CCF itself; b) times the formal RV error is a correct estimate of the error in an ideal case, that is, for a CCF without correlations among its points. Therefore, it is better to assume a more conservative estimate of the error. We increase the error on the bisector differences by a factor of , which corresponds to the assumption that there are only ten uncorrelated values along the CCF intervals used to estimate the BIS, and increase from to 2.5 the multiplicative factor applied to the RV error as given by the DRS. This factor is slightly larger than the generally adopted factor of 2.0 (cf. Sect. 4 of Boisse et al. 2009), yielding a conservative estimate of the BIS uncertainty. The maximum of those two estimates is assumed as our error on the BIS. An exact estimate of the number of independent bins along the CCF profile is difficult to obtain, so we follow Figueira et al. (2015) who estimate an oversampling factor of ~3 on the average (see below). In Sect. 4, we shall compare the errors on the BIS computed with the above procedure with the residuals of the regressions of the BIS time series to quantify the level of error overestimation introduced by our assumptions.
To compute the indicator ΔV of Nardetto et al. (2006) according to the definition in Sect. 5.1 of Figueira et al. (2013), a LevenbergMarquardt minimization algorithm is applied to fit a biGaussian and a Gaussian profiles to the normalized CCF, respectively. Only the portion of the CCF below 0.95 of the continuum is fitted to reduce the contribution of the Lorentzian wings of the photospheric lines to the best fit. The biGaussian fit has four free parameters: depth, central RV, FWHM, and an asymmetry measure. The Interactive Data Language (hereafter IDL) procedure mpfit.pro^{3} is used to compute the best fits. The indicator ΔV is defined as the difference in the central RV values obtained with the Gaussian and the biGaussian best fits, respectively. With the Gaussian best fit, we also compute the FWHM of the CCF and check that its value is within 0.1 percent of that provided by the DRS.
The errors on the parameters of the Gaussian and biGaussian best fits, from which ΔV and the FWHM of the CCF are derived, cannot be properly estimated from the covariance matrixes given by mpfit.pro because the shape of the CCF is not a Gaussian, which leads to systematic residuals when we perform a Gaussian best fit. Those systematic residuals are generally greater than the random and the correlated errors of the CCF and produce an overestimate of the errors on ΔV and the FWHM by a factor of five to seven when we use the standard method based on the covariance matrixes to evaluate their errors. Therefore, we decided to apply a different approach that takes into account both the effect of the random noise and the correlated noise.
The standard deviation σ_{A0(i)} of the CCF in the ith flux bin can be conservatively assumed to be , where A_{0}(i) is the value of the CCF in the bin. This happens because of the method applied to compute the CCF. Specifically, in the photondominated regime, the standard deviation of a given spectral element of the spectrum with N _{ph} photon counts is . When we compute the CCF value in a given bin, we add the photon counts in different spectral elements weighted according to the line mask adopted to compute the crosscorrelation. Since the mask weights are smaller than the unity, A _{0}(i) is smaller than the sum of the photoelectron counts in the spectral elements contributing to that bin of the CCF, implying that the signaltonoise ratio of the CCF is reduced with respect to that corresponding to the total photoelectron counts. Therefore, our prescription to compute σ_{A0(i)} is conservative.
To evaluate the effect of the correlated noise along the CCF, first we compute a smoothed version of the CCF itself, using a Savitzy–Golay filter of order four with ten points on the left and ten on the right of each point as explained in Ch. 14.9 of Press et al. (2002). The residuals between the CCF and its smoothed version are used to compute 100 realizations of the CCF with correlated noise by means of the prayerbead method (see Sect. 4.2 in Cowan et al. 2012). Computing ΔV and the FWHM for all of these 100 realizations, we can evaluate their standard deviations as produced by the correlated noise. In addition, we compute another 100 realizations by summing to the normalized CCF one hundred realizations of a Gaussian random noise of zero mean and standard deviation . Finally, we sum in quadrature the standard deviations of ΔV and the FWHM for the 100 realizations with correlated noise and the 100 realizations with random noise, thus obtaining their standard deviations that include the effects of both the correlated and the random noises.
The V_{asy} indicator of Figueira et al. (2013) incorporates the radial velocity RV(i) of each CCF element into its definition. Therefore, it changes if the RV scale is shifted by a constant amount or if the star is genuinely Doppler shifted due to the presence of a planet, thus giving rise to a misleading rejection of any true orbital motion of the star. Moreover, the value of V_{asy} is not invariant when the number of photoelectron counts along the spectrum varies, as is the case when a different exposure time is adopted or the star is observed at different air masses. These drawbacks were addressed by Figueira et al. (2015) who revised the definition of the indicator. Here, we propose an improved definition of V_{asy} to eliminate such spurious dependences and make the indicator nondimensional, that is,
where the weight at flux level i on the red or the blue wings of the CCF is defined as
where A_{0}(i) is the crosscorrelation value at flux level i and ∂A_{0}(i)/∂RV(i) the derivative of the CCF. The derivative is computed using a Savitzy–Golay filter of order four with ten points on the left and ten on the right of each given point. The mean weight for a given flux level is the arithmetic mean of the corresponding weights on the red and blue wings of the CCF. To avoid sampling too close to the continuum or to the bottom of the CCF, the summations in Eq. (1) are extended from five percent to 95 percent of the depth D of the CCF. The square of the weights in the denominator of Eq. (1) makes V_{asy(mod)} invariant for a multiplication of the number of the photoelectron counts by a constant factor.
To compute the error of V_{asy(mod)}, we consider both the effects of a Gaussian noise and of a correlated noise following the procedure previously described for ΔV and the FWHM. The final value of V_{asy} is the median over the 200 realizations of the CCF obtained with random and correlated noises. The square root of the sum of the variances of the 100 realizations with correlated noise and of the 100 with Gaussian noise is adopted as the standard deviation of the indicator, as in the case of ΔV and the FWHM. In light of the above considerations, this is a conservative estimate of the error on our indicator.
Figueira et al. (2015) consider the effect of the correlated noise in the CCF in the evaluation of their new V_{asy} indicator. Our V_{asy(mod)} is the same as their new indicator provided that, in their notation , that is, that the error on the RV as induced by the flux measurement in the bin i of the CCF is proportional to the value of the crosscorrelation in that bin. They recommend choosing the flux grid in such a way that it projects onto a grid of bins separated on the average by α_{corr} bins of the original CCF, where α_{corr} = (size of the bin in RV)/(step used in the CCF calculation) (see their Sect. 2.1 for details). In the case of the HARPS spectra, α_{corr} = 3.28. This implies that we oversample the CCF by the same factor when we compute our V_{asy(mod)} by considering all the flux levels between five and 95 percent of the continuum. We can correct for this effect by modifying our definition of the weight at flux level i as
Since the factor α_{corr} is constant and our definition of V_{asy(mod)} is invariant for any constant scale factor applied to the CCF value A_{0}(i), our value of the indicator is independent of α_{corr}. Therefore, we maintain our previous definition of V_{asy(mod)} and do not introduce any correction for the oversampling associated with the correlated nature of the CCF. We provide an IDL procedure to compute the asymmetry indicators and the FWHM of the CCF, which is described in Appendix A.
3.2. Correlations of the activity indicators with the RV variations
Our RV residual time series have no significant periodicities, therefore we cannot compare their periodograms with those of the activity indicators to investigate the correlation between RV and activity as for example in Santos et al. (2014), because our periodograms are dominated by random noise. Moreover, previous investigations showed that the correlations between the RV variations and our indicators generally are neither linear nor monotonic; for example, a plot of the RV vs. the BIS generally displays an eightshaped pattern (Boisse et al. 2011; Figueira et al. 2013). From this point of view, we are far from the ideal case of linearly correlated or even monotonically correlated indicators, thus Pearson or Spearman correlation coefficients are of very limited use.
Recently, Gaussian process (hereafter GP) regression has been proposed to model the RV variations produced by stellar activity (e.g. Haywood et al. 2014; Rajpaul et al. 2015). It is a nonparametric regression method that allows us to model complex time variations with both stochastic and deterministic components by parametrizing the covariance between pairs of data points. The covariance is specified by a set of hyperparameters that can be estimated from the dataset within a Bayesian framework (Roberts et al. 2012).
In the case of RV variations induced by stellar activity, the covariance is specified by a set of four hyperparameters associated with the amplitude of the correlation, the timescale for the growth and decay of active regions, the rotation period of the star, and the smoothing of the rotational modulation. These hyperparameters can be estimated either from RV time series with a suitable cadence and extension (e.g. LópezMorales et al. 2016) or using external datasets with high cadence, for example highprecision photometry (e.g. Malavolta et al. 2018).
In the present work, we explore an alternative approach to model the regression between the RV variations and the activity indicators. Specifically, we model the dependence of the RV on time and activity indicators by means of kernel regression (hereinafter KR), another nonparametric regression technique. Its foundations are given in, for example, Hardle (1992), Loader (1999), and Takeda et al. (2007), while previous applications to the analysis of astronomical data can be found in Wang et al. (2007), Marco et al. (2015), or AL Otaibi et al. (2016).
We apply a locally linear model to fit the RV at a time t_{k} by minimizing the objective function Z:
where N is the number of RV and simultaneous activity indicator observations, RV(t_{i}) the RV at time t_{i}, x_{i} ≡ x(t_{i}) a generic indicator at time t_{i}, β_{0} and β_{1} the coefficients with respect to which the objective function Z is minimized, and W the kernel given by
where h_{x} and h_{t} are the bandwidths of the kernel. Equation (4) is a linear regression of the RV vs. the activity proxy x where the datapoints are weighted according to the kernel in Eq. (5). This gives more weight to points closer in time to t_{k} and to the indicator x(t_{k}), thus accounting for the temporal variation of the correlation and the nonlinear dependence of the RV on the indicator itself. By the standard method to compute a weighted linear best fit (e.g. Press et al. 2002), we find the kernel regression estimator for the radial velocity at the time t_{k} as
where the elements of the matrix M are given by
where
with l = 0, 1, 2 and
The optimal values of the bandwidths h_{x} and h_{t} are obtained by the socalled leaveoneout method (Hardle 1992), that is, by minimizing the function
with respect to h_{x} and h_{t}, where the summation is made on the N_{r} datapoints that have a distance in x and t of at least 2h_{x} and 2h_{t} from the extreme points of the ranges in x and t, respectively, to avoid the effect of the worse quality of the fit at the border of the domain^{4}. We impose that h_{x} and h_{t} be smaller than 1/8 of the total range of the indicator x and of the time interval covered by the observations, respectively, to have a sufficient resolution of the kernel. Moreover, we impose that h_{x} and h_{t} be greater than 4 times the mean separation of the datapoints in x and t to avoid a too small number of datapoints effectively contributing to the regression at a given point. To perform the constrained minimization, we make use of the IDL procedure CONSTRAINED_MIN.
In addition to fitting the RV observations, KR can be used to interpolate the RV values between the times of the actual observations. To this purpose, first we linearly interpolate the values of x on an evenly sampled time grid {t_{j}} and then compute the kernel regression over the couples (t_{j}, x(t_{j})) to obtain a plot without time gaps. Confidence intervals for the estimator are computed by the method in Sect. 2.3.3 of Loader (1999).
The significance of the KR with respect to a given indicator x and the time t can be derived by a suitable application of the Fischer–Snedecor F statistics as discussed in Ch. 9 of Loader (1999). Specifically, we compute the ratioofvariance statistics F to compare the kernel regression model with the model assuming no dependence of the RV on the given indicator as
where is the mean of the RV measurements and is the effective number of degrees of freedom of the KR that is given by (cf. Loader 1999)
where Tr gives the trace of the argument matrix and Λ^{T} is the transpose of the matrix Λ that is defined as Λ ≡ I − M, with I being the identity matrix. The statistics F is distributed as the ratio of two χ^{2} variables with ν − 1 degrees of freedom in the numerator and N − ν − 1 degrees of freedom in the denominator, respectively. Its significance can be computed, for example, by means of the F_PDF function of IDL.
4. Results
The results of the application of the KR to the time series of the RV residuals obtained with the model including the keplerian motions of the known planets, a longterm linear trend, and a linear correlation between the RV and the chromospheric index are presented in Table A.1 at the end of the paper. A selection including only the cases where the significance of the KR is below 0.01 is presented in Table 2. We considered only the datapoints that had simultaneous measurements of the RV and of all the five indicators, that is , BIS, ≤V, V_{asy(mod)}, and FWHM, and applied the KR twice, first considering all the datapoints in the time series, then excluding the datapoints that deviated from the regression by more than three standard deviations of the residuals of the first regression. Tables 2 and A.1 list from left to right, the name of the star, the indicator used together with the time to compute the KR (cf. Eqs. (4) and (5)), the number N_{c} of datapoints after the above crossmatching and 3σ clipping, the standard deviation σ of these RV residuals before the application of the KR, the standard deviation σ_{KR} after the subtraction of the KR, the time bandwidth h_{t}, the indicator bandwidth h_{x}, the effective number of degrees of freedom of the KR, the Fischer–Snedecor function F that compares the regression model with that assuming no correlation between the RV residuals and the indicator and time, and the significance of the correlation, that is, the probability that a value of F as large as that observed can arise from random statistical fluctuations (cf. Sect. 3.2).
Results of the KR of the RV residuals of the keplerian model including a linear correlation between the RV and the stellar chromospheric index for our stars.
The results of the KR as applied to the RV residuals of the model including only the keplerian motions and a longterm linear trend in the RV are listed in Table A.2 at the end of the paper, which lists the same quantities as in Table 2. A selection of the results with α < 0.01 are reported in Table 3.
Results of the KR of the RV residuals of the keplerian model that does not include the linear correlation between the RV and the stellar chromospheric index for our stars.
The standard deviations of the residuals of the orbital best fits including or excluding the correlation between the RV and the chromospheric index are almost the same, except for HD 23596, HD 188015, and, to a lesser extent, HD 75898, for which they are smaller for the former model because of the additional free parameters. This confirms that the bias introduced by neglecting the dependence of the RV on activity is generally small for the stars in our sample^{5}. However, including the linear correlation makes a significant difference in those three cases. Stars HD 23596 and HD 188105 show a significant correlation between the RV variation and the activity indicators (α < 0.01; see Tables 2 and 3) suggesting the importance of including such a correlation in the keplerian fits. Thus we mainly refer to models including the correlation in the presentation of our results. Considering that we test five different correlations for 15 stars, that is a total of 75 possible correlations, we fix the significance threshold at α = 0.01 to have less than one expected spurious correlation arising from statistical fluctuations in our samples. Including the correlation in the keplerian model gives seven stars out of 15 with α < 0.01, while excluding that correlation gives six stars because the minimum α for HD 155358 increases from 0.0057 to 0.0247. Looking at Table 2, we see that the indicator with the greatest number of significant correlations is the BIS with four cases followed by the new indicator V_{asy(mod)} and the FWHM with two cases, while the chromospheric index and ΔV have only one. The significance of the KR of HD 75898 with respect to the V_{asy(mod)} and the time is 1.36 percent (see Table A.1). If we include this further case among the significant correlations, V_{asy(mod)} gives three positive cases. On the other hand, in the case of the keplerian model without the correlation, the BIS shows four significant correlations, the V_{asy(mod)} three, the FWHM two, and the ΔV and only one significant correlation (cf. Table 3).
In all the cases, the application of the KR reduces the amplitude of the RV residuals with their standard deviation ranging between 1 and 2 m s^{−1}, except for HD 75898, our most active star. Therefore, KR proved to be a suitable regression method to account for the effects of stellar activity on RV variations. In the cases of our significant correlations, the standard deviation of the KR residuals is always lower than 1.5 m s^{−1}, except for HD 188015, making a highly significant improvement in the modelling of the RV time series.
Two examples of the application of KR to our stars with well sampled time series are plotted in Figs. 1 and 2 for HD 106252 and HD 190228, respectively, both with 43 datapoints. The slope of the regression changes remarkably, sometimes discontinuously, at times when there are no datapoints to constrain the slope itself or when there are large variations in the RV residuals. Slope changes are more gradual in the time intervals when datapoints are more abundant and show variations with a lower amplitude. This is a consequence of the local linear fit performed on the datapoints and the need to match fits with different slopes. The performance of the regression is generally very good showing that it can capture the nonlinear and nonmonotonic correlations between the RV variations and the different activity indicators. We note that HD 106252 and HD 190228 are targets with a chromospheric index comparable with that of the Sun at the minimum of the elevenyear cycle.
Fig. 1. Top panel: radial velocity (RV) residuals of HD 106252 vs. time. The solid line is the kernel regression (KR) of the RV time series with respect to the time and the BIS; the ±σ interval of the regression as discussed in Sect. 3.2 is indicated (dotted lines). Bottom panel: RV residuals of the KR plotted in the top panel. Zero residuals are indicated by the dotted line. 
Fig. 2. Same as Fig. 1, but for HD 190228 with the KR performed with respect to the time and V_{asy(mod)}. 
Two examples of the application of the KR to stars with a relatively low number of datapoints are shown in Figs. 3 and 4 for HD 89307 with 19 points and HD 188015 with 21 points, respectively. The number of the degrees of freedom of the regression model is approximately half the number of the datapoints N_{c} as in the case of the stars with larger datasets indicating that the number of free model parameters N_{c} − ≈ 0.5 N_{c} in both the cases.
Fig. 3. Same as Fig. 1, but for HD 89307 with the KR performed with respect to the time and the V _{asy(mod)}. 
Fig. 4. Same as Fig. 1, but for HD 188015 with the KR performed with respect to the time and the FWHM. 
In principle, the bandwidths of our KR can be related to physical properties of the star. For example, h_{t} can be a measure of the evolutionary timescale of the pattern of active regions responsible for the RV variation, or h_{x} can be used to characterize a dependence of the RV perturbation on the level of activity. However, a word of caution is in order here because the sparseness of the sampling of our time series can limit the information on the active region lifetime and the dependence of the RV on the activity level that can be extracted from our datasets. Considering time series with more datapoints and less gaps such as those of HD 106252 or HD 190228, we estimate an average lifetime of the active regions from h_{t} between ~40 and ~80 days, which compares well with the typical lifetimes of solar faculae that dominate the RV variation in the case of the Sun (cf. Meunier et al. 2010; Lanza et al. 2016; Haywood et al. 2016).
In Sect. 3.1 we introduced a prescription to estimate the error ε on the BIS that was designed to overestimate the error itself. We can compare our BIS errors with the residuals of the regressions of the BIS time series of our stars to quantify the degree of overestimation. For each of our stars, we compute the residuals of the KR of the BIS vs. the time and the RV residuals, and consider the distribution of the ratio δ/ε for our whole sample of measurements. Ideally, such a distribution is expected to be a Gaussian with zero mean and unity standard deviation. Actually, we find that we must multiply ε by a correction factor γ = 0.44 to obtain a distribution with unity standard deviation as plotted in Fig. 5. We see that the fraction of measurements with very small deviations is significantly greater than expected in the case of a Gaussian distribution. Therefore, the Gaussian best fit plotted in the figure was computed by excluding the two central bins where virtually all the residuals with anomalously small deviations are concentrated. Those anomalous residuals amount to approximately ten percent of the total and are a consequence of the overfitting produced by the KR when the number of datapoints is low and their cadence is sparse. Nevertheless, they do not represent a major problem because their relative fraction is rather small. Alternatively, to avoid any problem with overfitting, we computed a simple linear regression of each of the BIS time series vs. the time only. The quality of those regressions is significantly worse than in the case of the KR, but we see no overabundance of residuals in the bins closer to zero. We find that the distribution of the ratio δ/ε can be well reproduced by a Gaussian with zero mean and a unity standard deviation when we adopt γ = 0.66 (see Fig. 6).
Fig. 5. Distribution of the ratio δ/ε of the KR residuals of the BIS timeseries of our stars vs. the time and the RV residuals with γ = 0.44. The solid line is a Gaussian best fit with unity standard deviation computed after excluding the two central bins from the fitting. 
Fig. 6. Distribution of the ratio δ/ε of the residuals of the linear regressions of the BIS timeseries of our stars vs. the time with γ = 0.66. The solid line is a Gaussian best fit to the distribution with unity standard deviation. 
We conclude that the procedure introduced in Sect. 3.1 overestimates the error on the BIS and that the average correction factor to be adopted to evaluate the true standard deviation is 0.44 < γ < 0.66. However, we do not introduce this correction factor in our procedure because we prefer to overestimate the errors on the BIS rather than underestimate them, given that γ is an average factor obtained from a sample including different stars. If required, in the case of a specific star having a homogenous and wellsampled dataset, the above statistical analysis can be repeated to correct a posteriori the BIS error obtained with our procedure.
5. Discussion and conclusions
We explored the performance of different asymmetry indicators of the CCF, its FWHM, and the chromospheric index to detect activityrelated RV variations in a sample of 15 slowly rotating (υ sin i ≲ 5 km s^{−1}), lowactivity () solarlike stars observed with HARPSN for two years. We gave prescriptions to compute CCF indicators and their errors starting from the CCFs provided by the DRS of HARPSN. After removing the RV modulations produced by known planets, we looked for periodic signals in the residuals finding no significant case for additional planets in our sample. Therefore, we computed a kernel regression (KR) of the residual RV with respect to the time and each of our activity indicators, including the newly defined V_{asy(mod)}, detecting at least one significant correlation (twosided pvalue < 0.01) in seven out of 15 stars. This gives an estimate of the global performance of the KR with our five indicators that was successful in ~47 ± 18 percent of the cases. The bisector inverse slope (BIS), the FWHM, and the new V_{asy(mod)} proved to be the most useful indicators because they provided a significant regression in four, two, and two cases out of 15, respectively.
In addition to the slow stellar rotation and low activity level, our relatively low performances could be a consequence of the conservative significance threshold (α < 0.01) chosen to minimize the occurrence of spurious correlations in our rather large sample. The performance of the CCF asymmetry indicators improves if we consider stars with an higher activity level or υ sin i such as HD 189733 or τ Bootis. With υ sin i = 2.97 ± 0.2 km s^{−1} and ranging from −4.524 and −4.458 (Pace 2013), HD 189733 shows a strong correlation between the BIS and the RV residuals with a falsealarm probability lower than 0.1 percent (Boisse et al. 2009). A similarly strong correlation has been detected in τ Boo that has υ sin i = 14.27 ± 0.06 km s^{−1} and ranging from −4.790 and −4.768 in the observations performed by Borsa et al. (2015). Our results demonstrate that KR can be fruitfully applied to our datasets that are characterized by a rather sparse sampling. Kernel regression offers advantages over simple linear or rank correlation analyses because it is not limited to assume a linear or monotone relationship between the RV variations and the activity indicators.
We note that KR uses two bandwidths that are estimated by minimizing the sum of the squared residuals of the regression (cf. Sect. 3.2). They account for the evolution timescale of the surface active regions and the nonlinear dependence of the RV perturbation on the activity indicators. In the implementation presented in this work, they are held fixed because we considered a time span of about two years, but it is simple to introduce a variation with the time to account for systematic changes along a stellar activity cycle with a typical duration of several years. Another advantage of the KR approach is the possibility of analytically estimating the significance of the regression by means of FisherSnedecor statistics.
A consistent comparison of KR with other techniques requires its inclusion in a global model of the RV variations taking into account both the keplerian motions and the effects of activity and allowing a Bayesian estimate of the most probable number of planets and of their orbital parameters (cf. Dumusque et al. 2017). This is outside the scope of the present paper, where we considered the evaluation of different asymmetry indicators of the CCF and adopted a simple approach to demonstrate the performance of KR with respect to different activity indicators considering a rather large sample of sunlike stars. The inclusion of KR into a complete RV model will be presented and discussed in future works dedicated to the modelling of the planetary systems of specific targets. The present results suggest that KR can account for a large part of the RV variations induced by stellar activity as shown by the remarkable reduction of the RV variability when we use KR with our suite of activity indicators. Our approach leads to a reduction of the standard deviations of the RV residuals that can exceed a factor of approximately two in the case of significant correlations (cf. Tables 2 and 3).
The best evaluation of the indicators discussed in the present work requires a highresolution and a stabilized spectrograph such as HARPSN. A spectral resolution allowing / of at least 10^{5} is recommended in the case of our slowly rotating stars as pointed out for example by Dumusque et al. (2014). With a lower resolution, the evaluation of the indicators is not optimal because the rotational broadening of the spectral lines cannot be resolved in stars with υ sin i ≲ 2−3 km s^{−1}. The longterm stability of the spectrograph is another fundamental issue because it warrants that the observed variation of the FWHM is not dominated by instrumental effects. Actually, we saw the remarkable variation of the FWHM produced by a slightly variable defocussing of HARPSN during the initial phase of its operation (see Benatti et al. 2017). On the other hand, the impact on the asymmetry indicators was found to be negligible in comparison with the errors associated with the photon shot noise at least in HARPSN. If its behaviour can be used as a guideline for similar highresolution spectrographs, this result suggests asymmetry indicators can be used instead of the FWHM when a high stability of the spectrograph is not warranted.
Available at https://github.com/LucaMalavolta/PyORBIT.
According to Dumusque et al. (2011b), ranges from −5.0 at minimum to −4.75 at maximum along the activity cycle of the Sun.
The justification of Eq. (10) can be found in, e.g., Loader (1999) or Sect. 5.2 of the notes by R. Tinshirani at http://www.stat.cmu.edu/~larry/=sml/nonpar.pdf.
Our stars are characterized by low levels of activity and rather large orbital RV modulations by their known planets, with orbital periods far from the typical rotation periods of sunlike stars. Specifically, the radial velocity semiamplitude K ranges between ~15 and ~140 m s^{−1}, while the orbital periods are between ~120 and ~4200 days. The activityinduced RV variations have most of their power at the rotation frequency and its first two or three harmonics (e.g. Boisse et al. 2011), which is far from the orbital frequencies and their harmonics for almost all of the known planets orbiting our stars. Therefore, the component of the activityinduced RV variation that can be absorbed by the orbital fit is small in most of our cases.
Acknowledgments
The authors wish to thank the anonymous referee whose valuable comments on a previous version of this work allowed them to improve the methodology of data analysis and the results. The present work is based on observations made with the highresolution spectrograph HARPSN at the Italian Telescopio Nazionale Galileo (TNG). Data have been made available through the INAFTNG IA2 national archive at INAFOsservatorio astronomico di Trieste whose availability is grateful acknowledged. Simulations obtained with SOAP 2.0 through its web interface (Dumusque et al. 2014) are gratefully acknowledged. The research leading to these results received funding from the European Union Seventh Framework Programme (FP7/20072013) under grant agreement number 313014 (ETAEARTH). Support from the Progetti Premiali scheme (Premiale WOW) of the Italian national Ministry of Education, University, and Research is also acknowledged. PF acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contracts of reference IF/01037/2013/CP1191/CT0001, respectively, and POPH/FSE (EC) by FEDER funding through the programme "Programa Operacional de Factores de Competitividade  COMPETE”. PF further acknowledges support from FCT in the form of an exploratory project of reference IF/01037/2013/CP1191/CT0001. The authors are grateful to Dr. R. Zanmar Sanchez for his comments on the present work.
References
 AL Otaibi, S., Tiňo, P., CuevasTello, J. C., Mandel, I., & Raychaudhury, S. 2016, MNRAS, 459, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benatti, S., Desidera, S., Damasso, M., et al. 2017, A&A, 599, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Biazzo, K., Gratton, R., Desidera, S., et al. 2015, A&A, 583, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boisse, I., Moutou, C., VidalMadjar, A., et al. 2009, A&A, 495, 959 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boisse, I., Bouchy, F., Hébrard, G., et al. 2011, A&A, 528, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borsa, F., Scandariato, G., Rainer, M., et al. 2015, A&A, 578, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2012, Proc. SPIE, 8446, 84461V [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2014, Proc. SPIE, 9147, 91478C [Google Scholar]
 Covino, E., Esposito, M., Barbieri, M., et al. 2013, A&A, 554, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cowan, N. B., Machalek, P., Croll, B., et al. 2012, ApJ, 747, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Damasso, M., Biazzo, K., Bonomo, A. S., et al. 2015a, A&A, 575, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Damasso, M., Esposito, M., Nascimbeni, V., et al. 2015b, A&A, 581, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desidera, S., Sozzetti, A., Bonomo, A. S., et al. 2013, A&A, 554, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desidera, S., Bonomo, A. S., Claudi, R. U., et al. 2014, A&A, 567, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desort, M., Lagrange, A.M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Díaz, R. F., Ségransan, D., Udry, S., et al. 2016, A&A, 585, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Donati, J.F., Hébrard, E., Hussain, G., et al. 2014, MNRAS, 444, 3220 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X. 2014, ApJ, 796, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011a, A&A, 525, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dumusque, X., Santos, N. C., Udry, S., Lovis, C., & Bonfils, X. 2011b, A&A, 527, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Glenday, A., Phillips, D. F., et al. 2015, ApJ, 814, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Borsa, F., Damasso, M., et al. 2017, A&A, 598, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Figueira, P., Santos, N. C., Pepe, F., Lovis, C., & Nardetto, N. 2013, A&A, 557, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Figueira, P., Santos, N. C., Pepe, F., Lovis, C., & Nardetto, N. 2015, A&A, 582, C2 [CrossRef] [EDP Sciences] [Google Scholar]
 Fischer, D. A., AngladaEscude, G., Arriagada, P., et al. 2016, PASP, 128, 066001 [NASA ADS] [CrossRef] [Google Scholar]
 Han, E., Wang, S. X., Wright, J. T., et al. 2014, PASP, 126, 827 [NASA ADS] [CrossRef] [Google Scholar]
 Hardle, W. 1992, Applied Nonparametric Regression, (Cambridge: Cambridge University Press) available at: http://ftsipil.unila.ac.id/dbooks/applied%20nonparametric%20regression.pdf [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [NASA ADS] [CrossRef] [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Unruh, Y. C., et al. 2016, MNRAS, 457, 3637 [NASA ADS] [CrossRef] [Google Scholar]
 Hébrard, É. M., Donati, J.F., Delfosse, X., et al. 2014, MNRAS, 443, 2599 [NASA ADS] [CrossRef] [Google Scholar]
 Lagrange, A.M., Desort, M., & Meunier, N. 2010, A&A, 512, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Boisse, I., Bouchy, F., Bonomo, A. S., & Moutou, C. 2011, A&A, 533, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Molaro, P., Monaco, L., & Haywood, R. D. 2016, A&A, 587, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Loader, C. 1999, Local Regression and Likelihood (New York: SpringerVerlag) [Google Scholar]
 LópezMorales, M., Haywood, R. D., Coughlin, J. L., et al. 2016, AJ, 152, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Lovis, C., Dumusque, X., Santos, N. C., et al. 2011, ArXiv eprints, [arXiv:1107.5325] [Google Scholar]
 Malavolta, L., Nascimbeni, V., Piotto, G., et al. 2016, A&A, 588, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malavolta, L., Mayo, A. W., Louden, T., et al. 2018, AJ, 155, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Marco, F. J., Martínez, M. J., & López, J. A. 2015, AJ, 149, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Bonfils, X., Forveille, T., et al. 2009, A&A, 507, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meunier, N., Desort, M., & Lagrange, A.M. 2010, A&A, 512, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nardetto, N., Mourard, D., Kervella, P., et al. 2006, A&A, 453, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pace, G. 2013, A&A, 551, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical recipes in C++ : the art of scientific computing, ed. H. William (Cambridge: Cambridge University Press) [Google Scholar]
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Queloz, D., Bouchy, F., Moutou, C., et al. 2009, A&A, 506, 303 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269 [NASA ADS] [CrossRef] [Google Scholar]
 Roberts, S., Osborne, M., Ebden, M., et al. 2012, Phil. Trans. R. Soc. A, 371, 20110550 [Google Scholar]
 Robertson, P., Endl, M., Cochran, W. D., et al. 2012, ApJ, 749, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Santos, N. C., Sousa, S. G., Mortier, A., et al. 2013, A&A, 556, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Santos, N. C., Mortier, A., Faria, J. P., et al. 2014, A&A, 566, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Takeda, H., Farsiu, S., & Milanfar, P. 2007, IEEE Transactions on Image Proc., 16, 349 [NASA ADS] [CrossRef] [Google Scholar]
 Udry, S., Mayor, M., Benz, W., et al. 2006, A&A, 447, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Unruh, Y. C., Solanki, S. K., & Fligge, M. 1999, A&A, 345, 635 [NASA ADS] [Google Scholar]
 Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Valenti, J. A., Fischer, D., Marcy, G. W., et al. 2009, ApJ, 702, 989 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, D., Zhang, Y. X., Liu, C., & Zhao, Y. H. 2007, MNRAS, 382, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., Kürster, M., & Endl, M. 2009, A&A, 505, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: An IDL procedure to compute the asymmetry indicators
We provide a procedure written in IDL 8.4 to compute the CCF asymmetry indicators and FWHM according to the methods described in Sect. 3. It provides also the old indicator V _{asy} as defined by Figueira et al. (2013) for reference to previous work, although it is not recommended due to its tendency to show spurious correlations with the RV variations.
To compile our procedure, installation of the IDL Astronomy Library is required together with the IDL procedure mpfit.pro (we used version 1.82) that can be downloaded from http://cow.physics.wisc.edu/~craigm/idl/fitting.html with its documentation. You will also need readcol.pro to read ASCII input files and readfits.pro that we included into our file for simplicity, although they are generally found in the IDL Astronomy Library.
The input data to our procedure consist of a set of fits files produced by the HARPS data reduction software (DRS), specifically those containing the CCF and the bisector profiles, the names of which are of the kind HARPN.YYYYMMDDTHHMMSS.SSS_ccf_MASK_A.fit and HARPN.YYYYMMDDTHHMMSS.SSS_bis_MASK_A.fit, respectively. The naming convention for indicating the reduced data files of the DRS is introduced in the HARPSN User Manual (Sect. 8) and in the DRS User Manual that are accessible through the web page: http://www.tng.iac.es/instruments/harps/. We recommend reading those manuals before using our IDL procedure. The bisector files are used only for the purpose of comparing our bisector profiles with those provided by the DRS.
The output of our procedure is an ASCII file giving for each of the input ✻_ccf_✻ and ✻_bis_✻ files the asymmetry indicators, the FWHM, and the contrast (central depth relative to the continuum) of the CCF together with their uncertainties. The correction factor γ for the BIS error, introduced in the final part of Sect. 4, is not included. Our procedure provides the standard deviations of ΔV and FWHM computed with the method described in Sect. 3 together with those obtained from the covariance matrix and the best fit residuals given by mpfit.pro, although the latter are not recommended for a proper evaluation of the uncertainties because they overestimate the errors (cf. Sect. 3). Details on the output can be found in the header comment lines of the procedure itself.
It is possible to use the procedure in an interactive way that allows the user to see screen plots (e.g. on a X11 terminal) of the CCF and its bisector as well as of the fitting functions used to evaluate the indicators. We recommend the analysis of a sequence of CCF profiles first interactively to see the actual contents of the dataset and check the proper fitting of the CCF profiles. After that stage, it is possible to run the procedure automatically to speed up the analysis of the dataset. When some best fits appear to be not completely adequate because of deviations exceeding a typical threshold of 0.1 percent, the procedure warns the user, asking whether to continue or not. In the case of a lack of sufficient information or bad fits, it stops printing a message. We enclose a set of ✻_ccf_✻ and ✻_bis_✻ files together with the corresponding output file to test the operation of the procedure.
A tarfile including the IDL macro and auxiliary files to compile and test its operation can be downloaded from the DropBox: https://www.dropbox.com/sh/74z4hb2wksqg9wg/AABak08BU3EjXFKr9XnSVYPva?dl=0 or from INAF GitLab at https://www.ict.inaf.it/gitlab/antonino.lanza/HARPSN_spectral_line_profile_indicators.git
Results of the KR of the RV residuals of the keplerian model including a linear correlation between the RV and the stellar chromospheric index for our stars.
Results of the KR of the RV residuals of the keplerian model that does not include the linear correlation between the RV and the stellar chromospheric index for our stars.
All Tables
Results of the KR of the RV residuals of the keplerian model including a linear correlation between the RV and the stellar chromospheric index for our stars.
Results of the KR of the RV residuals of the keplerian model that does not include the linear correlation between the RV and the stellar chromospheric index for our stars.
Results of the KR of the RV residuals of the keplerian model including a linear correlation between the RV and the stellar chromospheric index for our stars.
Results of the KR of the RV residuals of the keplerian model that does not include the linear correlation between the RV and the stellar chromospheric index for our stars.
All Figures
Fig. 1. Top panel: radial velocity (RV) residuals of HD 106252 vs. time. The solid line is the kernel regression (KR) of the RV time series with respect to the time and the BIS; the ±σ interval of the regression as discussed in Sect. 3.2 is indicated (dotted lines). Bottom panel: RV residuals of the KR plotted in the top panel. Zero residuals are indicated by the dotted line. 

In the text 
Fig. 2. Same as Fig. 1, but for HD 190228 with the KR performed with respect to the time and V_{asy(mod)}. 

In the text 
Fig. 3. Same as Fig. 1, but for HD 89307 with the KR performed with respect to the time and the V _{asy(mod)}. 

In the text 
Fig. 4. Same as Fig. 1, but for HD 188015 with the KR performed with respect to the time and the FWHM. 

In the text 
Fig. 5. Distribution of the ratio δ/ε of the KR residuals of the BIS timeseries of our stars vs. the time and the RV residuals with γ = 0.44. The solid line is a Gaussian best fit with unity standard deviation computed after excluding the two central bins from the fitting. 

In the text 
Fig. 6. Distribution of the ratio δ/ε of the residuals of the linear regressions of the BIS timeseries of our stars vs. the time with γ = 0.66. The solid line is a Gaussian best fit to the distribution with unity standard deviation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.