Separating planetary reflex Doppler shifts from stellar variability in the wavelength domain

ABSTRACT

Stellar magnetic activity produces time-varying distortions in the photospheric line profiles of solar-type stars. These lead to systematic errors in high-precision radial-velocity measurements, which limit efforts to discover and measure the masses of low-mass exoplanets with orbital periods of more than a few tens of days. We present a new data-driven method for separating Doppler shifts of dynamical origin from apparent velocity variations arising from variability-induced changes in the stellar spectrum. We show that the autocorrelation function (ACF) of the cross-correlation function used to measure radial velocities is effectively invariant to translation. By projecting the radial velocities on to a subspace labelled by the observation identifiers and spanned by the amplitude coefficients of the ACF’s principal components, we can isolate and subtract velocity perturbations caused by stellar magnetic activity. We test the method on a 5-yr time sequence of 853 daily 15-min observations of the solar spectrum from the HARPS-N instrument and solar-telescope feed on the 3.58-m Telescopio Nazionale Galileo. After removal of the activity signals, the heliocentric solar velocity residuals are found to be Gaussian and nearly uncorrelated. We inject synthetic low-mass planet signals with amplitude K = 40 cm s⁻¹ into the solar observations at a wide range of orbital periods. Projection into the orthogonal complement of the ACF subspace isolates these signals effectively from solar activity signals. Their semi-amplitudes are recovered with a precision of ∼ 6.6 cm s⁻¹, opening the door to Doppler detection and characterization of terrestrial-mass planets around well-observed, bright main-sequence stars across a wide range of orbital periods.

1 INTRODUCTION

For decades, Doppler spectroscopy has been one of the most productive methods to discover and characterize exoplanets. Improvements in the precision, wavelength calibration, and stability of high-resolution échelle spectrographs have allowed exoplanet surveys to probe planets with radial velocity (RV) semi-amplitudes of just ∼1 m s⁻¹. New generations of spectrographs such as CARMENES (Quirrenbach et al. 2014), ESPRESSO (Mégevand et al. 2014), EXPRES (Jurgenson et al. 2016), HARPS-3 (Thompson et al. 2016), HPF (Ninan et al. 2018), and NEID (Schwab et al. 2016) are being designed and commissioned with improved resolution, spectral coverage, wavelength calibration, and stabilization systems (Wright & Robertson 2017). Recently, ESPRESSO has achieved 30 cm s⁻¹ precision per RV observation on Proxima Cen (Suárez Mascareño et al. 2020), and EXPRESS has achieved 58 cm s⁻¹ precision per RV observation on HD 3651 (Brewer et al. 2020).

Even with present instruments, the ability of spectroscopic surveys to detect and characterize low-mass planets is often limited by stellar variability and the stability of the wavelength calibration, rather than photon noise or instrumental errors (e.g. Saar & Donahue 1997; Queloz et al. 2001; Haywood et al. 2014). The purpose of this study is to devise a practical new approach to measuring stellar RVs in a way that mitigates the errors due to line-shape changes caused by stellar variability. To achieve this, we make use of the fact that changes in the shape of spectral lines may influence the apparent RV, but changes in the range rate (the first derivative with respect to time of the distance from the star’s centre to the Solar-system barycentre) induce only a shift and do not affect the line shape or depth. Related approaches exploiting profile-shape changes of even and odd character to disentangle shifts from activity have been published recently by Zhao & Tinney (2020) and Holzer et al. (2020) while de Beurs et al. (2020) have employed a neural-network machine learning to relate activity-related RV shifts to cross-correlation function (CCF) profile-shape changes in the same solar data set examined here.

One feature of the method we present is that it makes use of existing data products, i.e. the CCFs between the observed spectra and a digital mask (see Section 2), and the RVs derived from them. The CCF C(v, t) is a function of both barycentric velocity v and time t. The temporal variability of the CCF includes both Doppler shifts and line shape changes. The derivatives of the CCF with respect to velocity [e.g. C′(v, t), C″(v, t)] are also functions of velocity and time.

Our long-term goal is to improve the detection sensitivity and robustness of spectroscopic planet searches. As a specific objective towards that goal, we aim to devise efficient, shift-invariant metrics that can contribute to characterizing the detailed line-profile shape at each epoch. We describe the building blocks for our new method in Section 2. Then, in Section 3, we propose a new algorithm to compute such metrics, based on a novel combination of existing data products and techniques. We apply these metrics as predictors for the contribution of stellar variability to the apparent Doppler shift and infer a cleaned velocity time series. In Section 4, we verify and validate the method based on injection and recovery tests using solar observations. While our data-driven method makes no assumptions about the physical origin of the stellar variability, the tests in this paper focus on magnetic activity of the Sun. Finally, we discuss the implications of our work for future spectroscopic planet surveys in Section 5.

2 CROSS-CORRELATION FUNCTION AND ITS AUTOCORRELATION

The HARPS, HARPS-N, and ESPRESSO data reduction systems (DRS) derive RVs from stellar échelle spectra by computing the CCF of the spectrum with a digital line mask matched approximately to the spectral type of the target star (Baranne et al. 1996; Pepe et al. 2002).

The data presented in this paper were re-analysed with the ESPRESSO data-reduction pipeline (Pepe et al. 2021). The reduction procedure differs from that used for previous analyses of HARPS-N data (cf. Collier Cameron et al. 2019) in several important respects. These are described in detail by Dumusque et al. (2020). The wavelength scale is derived using a new line list tailored for the primary HARPS-N ThAr calibration lamp. A single master wavelength calibration is used for all observations. The wavelength of each pixel in the extracted spectrum is corrected for drift relative to the master calibration using the daily wavelength calibrations with the primary ThAr lamp, and the simultaneous reference source (either a secondary ThAr lamp or a Fabry–Perot). Prior to cross-correlation with a synthetic mask, the pixel wavelengths are transformed to the reference frame of the Solar system barycentre along the line of sight to the target. The CCFs are estimated on a common grid of pixels at uniformly spaced intervals of h = 0.82 km s⁻¹ in velocity space, using a blaze-corrected, inverse-variance weighted cross-correlation with a mask of line wavelengths and weights appropriate to the target spectral type. The CCF sampling interval matches the velocity increment per physical CCD pixel in the instrument. The velocity scale of the CCF is in the reference frame of the Solar system barycentre, with the drift correction applied.

Being the cross-correlation of a stellar absorption spectrum with a positive line mask, the CCF computed by the DRS resembles a single stellar absorption line, with a pseudo-continuum level which is accurately and consistently normalized to unity. The resulting CCF profile is fitted with one minus a Gaussian function described by three parameters: central velocity (v, relative to the mask), full width at half-maximum depth (FWHM) and central line depth as a fraction of the pseudo-continuum. Similarly, the bisector inverse slope (BIS) of the profile (Queloz et al. 2001) is recorded as a measure of profile asymmetry.

Stellar activity compromises the fidelity of this method of RV measurement. The contrast between bright convective upflows in photospheric granules and cooler downflows in intergranular lanes imposes an inherent asymmetry on the line profile (Dravins, Lindegren & Nordlund 1981). Magnetic activity, the finite lifetime of the granulation pattern and P-mode oscillations all cause the already non-Gaussian shape of the observed CCF to vary with time. As the star rotates, dark star-spots produce line-absorption deficits which migrate across the profile from blue to red, introducing time-varying amounts of skew and kurtosis into the spectral-line shapes (Saar & Donahue 1997; Dumusque, Boisse & Santos 2014). These magnetic regions alter the local convective velocity and line-profile asymmetry, as well as the local brightness weighting of the stellar rotation profile (Meunier, Lagrange & Desort 2010). In faculae-dominated stars like the Sun, magnetic suppression of granular convection in faculae causes even stronger time-varying profile asymmetries than sunspots, combining rotational Doppler shifts with foreshortening-dependent changes in the radial-tangential velocity field (Meunier et al. 2010; Cegla et al. 2019).

Several previous studies (e.g. Aigrain, Pont & Zucker 2012; Dumusque et al. 2012; Rajpaul et al. 2015) have explored whether the estimated velocities could be improved by decorrelating with other measurements, such as the FWHM or BIS. While the FWHM and BIS contain information about the CCF shape, they are not sufficient to describe the detailed changes in the CCF. In order to make use of all information contained in the spectra, some studies have suggested analysing stellar variability by applying a principal-component analysis (PCA) to the observed spectroscopic time series (e.g. Davis et al. 2017; Jones et al. 2017). These methods are promising on simulated data sets, but applying this approach to actual observations is challenging due to details of the spectrograph and its calibration process. In this paper, we apply the PCA approach to the CCF instead of the raw spectrum. This approach leverages extensive investments in developing a robust pipeline to measure the CCF. Of course, analysing only the CCF does reduce the total information content of the spectrum. We offer suggestions for how the method could be generalized to extract more information in Section 5.

To illustrate the impact of stellar variability on the CCF and to provide a test data set, we created a sequence of 886 daily solar CCFs from the HARPS-N solar telescope (Cosentino et al. 2014; Dumusque et al. 2015; Phillips et al. 2016), spanning the period from 2015 July 29 to 2020 March 6. The resulting distortions of the CCF profile, described in detail by Collier Cameron et al. (2019), are seen most clearly in the residuals |$R(v_i,t_j)= \mathbf {C}(v_i,t_j)- \left\langle \mathbf {C}(v_i)\right\rangle $| obtained by subtracting the time-averaged profile of the entire 5-yr sequence from each CCF in the sequence (Fig. 1, left-hand panels).

Full details of the HARPS-N solar observations are given by Collier Cameron et al. (2019), who used a Gaussian mixture model to assign a probability to each observation that it is unaffected by uneven transparency across the solar disc and corrected the velocities for differential extinction. Here we select only those observations with (1) probabilities greater than 99 per cent of being good and (2) with velocity corrections for differential extinction less than 10 cm s⁻¹. In summer these conditions are satisfied for up to 4 h d⁻¹, and in winter for up to 2 h d⁻¹.

Our 5-min exposure time is dictated by the need to average out solar p-mode oscillations. Light from the Sun is gathered by a 76 mm objective lens of 200 mm focal length, scrambled in an integrating sphere and fed into the spectrograph via an optical-fibre feed to the calibration unit and a neutral-density filter which attenuates the throughput by a further factor 15 (Phillips et al. 2016). The overall throughput is comparable to night-time HARPS-N exposures for a star of magnitude 5.5, for which we use the same exposure time in good seeing without saturating the detector. This gives signal-to-noise ratio (SNR) ≃ 350 or so in the continuum in échelle order 60, which translates to SNR ≃ 5000 in the CCF. For stellar RV observations we use 15-min blocks of contiguous exposures to mitigate the effects of p-mode oscillations. Within the windows that satisfy our selection criteria each day, we select at random a set of three contiguous CCFs to mimic an RV observation of a V = 5.5 star, and form a weighted average CCF using the square of the mean SNR of the CCF as the weighting factor. Thus, we anticipate that our test data set is dominated by solar magnetic activity, granulation and/or instrumental issues. We encourage future studies to investigate how well the algorithm can mitigate spectral variability on shorter time-scales.

2.1 Translation to the heliocentric frame

The HARPS-N DRS was designed primarily for stellar RV measurement, computing the CCF in the reference frame of the Solar-system barycentre in the direction of the target, as described above. As a result, the instantaneous CCF is Doppler-shifted by the component of the Sun’s barycentric motion in the observer’s direction, as is apparent from the upper-left panel of Fig. 1. To convert CCFs derived from solar spectra to the heliocentric reference frame, the CCF profiles must be shifted by the line-of-sight component ϵ of the Sun’s reflex motion about the barycentre.

The barycentric to heliocentric velocity corrections were computed using the JPL horizons software of Giorgini et al. (1996).

The derivatives are calculated numerically using equations (A1) and (A2) (see Appendix  A). The differences between neighbouring CCF values are substantially less than unity. The barycentric to heliocentric velocity correction is never greater than ±14.7 m s⁻¹, which is much less than the h = 820 m s⁻¹ sampling interval between neighbouring CCF elements. The truncation error in equation (1) is therefore significantly less than ϵ³/12h³ ≃ 2.6 × 10⁻⁸.

We validated the fidelity of the shift by calculating heliocentric velocities from the shifted profiles using the methodology of Appendix  B. We computed barycentric velocities from the original un-shifted CCFs by the same method, then applied the barycentric to heliocentric velocity correction. The RMS scatter of the difference between the two resulting sets of heliocentric velocities was 0.008 m s⁻¹. The RMS scatter in the difference between the heliocentric velocities calculated from the shifted CCFs and the DRS velocities transformed to the heliocentric frame was 0.044 m s⁻¹.

The resulting CCF time series, shown at lower left in Fig. 1, is effectively that of a star with no planets. This image shows that the form of the solar CCF is far from static, with dramatic changes in profile shape taking place on all time-scales from days to years.

2.2 Autocorrelation of CCF and shift invariance

The right-hand panels of Fig. 1 show that, despite the strong differences between the residual CCFs in the barycentric and heliocentric frames, their ACFs are very similar. The similarity is, however, only approximate. The autocorrelation domain is not infinite, and the circular shift method employed is vulnerable to edge effects if the CCFs are strongly shifted. In Fig. 2 we compare the column variances of the barycentric ACF time series with the column variances of the residuals obtained by subtracting the barycentric from the heliocentric ACF time series. We find that the temporal variance of the residual ACF is between 2.5 and 4.5 orders of magnitude smaller than the temporal variance of either the barycentric or heliocentric ACF, at every point in the profile. We conclude that for the purposes of this study, the ACF is effectively invariant to the solar reflex motion around the Solar-system barycentre.

2.3 Principal-component analysis of stellar variability

The columns of |$\mathbf {U}_{C,k}(t_j)$| and |$\mathbf {U}_{A,k}(t_j)$| define an orthonormal basis in the time domain. Each column comprises the coefficients (also known as scores) that define the temporal behaviour of the corresponding eigenvector. Fig. 3 (bottom) shows the scores of the leading 6 eigenvectors for all the individual observations in the time series ensemble, plotted against barycentric Julian date. The ACF is calculated in such a way that it is an even function, so its eigenvectors are also even functions. Those of the CCF display a mix of even and odd character. Nonetheless, there are strong similarities in the temporal behaviours of their scores.

The scores of the first principal component of both the CCF and the ACF, plotted in blue at the top of all four panels, show a secular upward trend with a superposed signal of higher frequency. The form of the trend, and the shape of the corresponding CCF eigenvector, indicates that this mode of variation affects both the depth and asymmetry of the line profile. It bears a strong resemblance to the variability of the CCF area (i.e. the product of the FWHM and central line depth) noted by Collier Cameron et al. (2019). These authors attributed the trend in CCF area to a secular decline in solar network flux and the faster variations to passages of active-region faculae across the solar disc. Thus, one sees that the ACF is able to recover a very similar series of scores in a way that is insensitive to line shifts. We will exploit this property for separating true Doppler shifts from stellar variability in Section 3.

The time variations of the scores of the second principal component of the ACF (orange traces, second from top in right-hand panels of Fig. 3) and the fourth principal component of the CCF (red traces, fourth from top in left-hand panels) are also similar, though the CCF version appears noisier. Collier Cameron et al. (2019) noted the same pattern of variability in the FWHM of the Gaussian profile fitted to the CCF by the HARPS-N DRS, arising from seasonal changes in the apparent solar rotational broadening. The Earth’s orbital eccentricity gives rise to an annual modulation in its orbital angular velocity, and hence the apparent solar rotation rate. The six-month oscillation in the obliquity of the solar rotation axis to the Earth’s orbital plane also affects the rotational broadening. The Bayesian Generalized Lomb–Scargle (BGLS; Mortier et al. 2015) periodogram of the second principal component of the ACF (Fig. 4) shows both periods clearly. The corresponding eigenvector for the CCF resembles the second derivative of the line profile, as expected for the CCF changing in width.

The third principal component of the ACF (green traces, third from top in right-hand panels of Fig. 3 and the CCF (green, third from top in left-hand panels) also resemble each other. They show apparently stochastic discontinuities followed by quasi-exponential decay with a time constant of order a few tens of days. These discontinuities are of instrumental origin. There is a very slow leak in the continuous-flow cryostat of the HARPS-N CCD. The cryostat has to be warmed up approximately twice per year to drive off the water that starts to obstruct the flow of liquid nitrogen. It has been observed that these warm-ups cause a sudden change in the asymmetry of the PSF, which takes a few weeks to decay (Dumusque et al. 2020). During the period of these observations, warm-ups were carried out at JD (24)57161.5, 57308.5, 57478.5, 57687.5, 57854.5, 58071.5, 58231.5, 58412.5, 58554.5, 58684.5, and 58839.5. These clearly coincide with the discontinuities in the third principal component of the ACF.

Examination of the BGLS periodograms (Fig. 4) of the leading principal components of the CCF reveals power at half the solar rotation period in the first, second, third, and fifth principal components. The ACF shows power at this period in the first, fourth, and fifth components. These components probably track profile-shape changes caused by sunspot groups and faculae traversing the visible solar hemisphere. There is surprisingly little power at the solar rotation period in the principal components of either the CCF or the ACF.

The CCF shows power at 6 months and/or 1 yr in its second and fourth components; the ACF shows power at these periods in the second and sixth components. Since the heliocentric time series by definition contains no solar reflex motion, these periodic shifts must also be associated with CCF profile-shape variability arising from Earth’s orbital motion.

The second principal component of the CCF and the third component of the ACF in Fig. 4, which we have identified with cryostat warm-ups, shows no power at the solar rotation period or its harmonics in either the CCF or the ACF, but we see significant structure on time-scales upwards of 100 d. This indicates that the changes in profile shape caused by cryostat warm-ups are different in character from any form of rotationally modulated solar activity.

Overall, we see that in the CCF, a simple profile shift would have non-zero projection on to multiple eigenvectors, including those that primarily represent broadening, skew and kurtosis. The eigenvectors of the corresponding components of the CCF have a mix of even and odd characteristics, and their odd parts should therefore affect the measured RVs. However, most of the time variations of the CCF appear broadly similar to those of the shift-invariant profile-shape changes probed by the ACF. This raises the possibility that the ACF can be used to deduce the contribution of profile shape changes to the measured RVs. We conclude that, at least for the heliocentric solar time series, principal-component analysis of the ACF could provide an effective means of separating the effects of dynamical shift from those of stellar and instrumental profile variability.

3 THE RV RESPONSE TO ACF TIME-DOMAIN VARIATIONS

To achieve this, we treat the time series of scores for each principal component of the ACF as the coefficients of a set of unknown eigenvectors representing orthogonal modes of variability in the shape of the CCF. These unknown eigenvectors of the CCF will affect the measured RV to a greater or lesser degree depending on whether they are pre-dominantly of even or odd character.

3.1 Projection into the ACF time-domain subspace

The set of RV observations has m elements, and can be thought of as a vector |$\mathbf {v}_{\rm obs}$| belonging to an m-dimensional space S. The l orthonormal columns of |$\mathbf {U}_A$| have the same dimension as S, and define an l-dimensional subspace U⊂S, centred on the origin of S.

From here on, we refer to |$\mathbf {v}_{\rm obs}$| as the ‘measured’ or ‘observed’ velocities’; |$\mathbf {v}_\parallel$| as the ‘model’ or ‘shape-driven’ velocities; and |$\mathbf {v}_\perp$| as the ‘shift-driven’ velocities. Fig. 5 shows that the shape-driven velocity perturbations, |$\mathbf {v}_\parallel$|⁠, are strongly correlated with the observed velocities, reproducing faithfully the long-term and short-term fluctuations dominated by stellar activity.

Given a set of measured RVs (⁠|$\mathbf {v}_{\rm obs}$|⁠) and the corresponding array of CCFs from which they were derived, equations (3), (5), (7), and (8) constitute a simple linear projection method for deriving the shape-driven perturbations to the RV (⁠|$\mathbf {v}_\parallel$|⁠), so as to provide a substantially cleaned set of shift-driven RVs (⁠|$\mathbf {v}_\perp$|⁠).

3.2 Outlier clipping

The shape of the CCF is sensitive to more than just solar activity. Changes in spectrograph focus can affect the FWHM of the CCF, while cryostat warm-ups perturb the skewness of the profile. Noisy CCFs, saturated exposures, or undetected cloud obscuration of part of the solar disc, can also cause temporary profile distortions which may not correlate with any of the highest-variance principal components.

Such anomalous observations may indeed generate their own basis functions when SVD is applied to the ACF time series. Their coefficients are normally close to zero, except when an anomaly occurs. They then appear as outliers in the corresponding columns of |$\mathbf {U}_A$|⁠. Such points can be masked as bad (0) if their absolute deviations lie further from the median value of the column than a specified number of median absolute deviations (MAD), and good (1) otherwise. If even one of the coefficients for an observation is an extreme outlier, it is likely that the entire observation is contaminated. We therefore create a one-dimensional rejection mask in the time domain from the product of the column masks. For the solar data we found that clipping at 6 times the MAD within each column of |$\mathbf {U}_A$| provided a stable set of basis vectors at the cost of reducing the total number of usable days of observation from 886 to 853.

This clipping procedure ensures a clean set of basis vectors, but does not detect outliers caused by unwanted velocity shifts, such as might be caused by an anomalous drift measurement. If present, these must be identified and clipped separately.

3.3 Rank reduction and column re-ordering

Following outlier clipping and masking, the singular-value decomposition of both the CCF and the ACF is re-computed from the surviving observations. It should be noted that all figures in this paper from Fig. 1 onward are based on the masked data set only.

The subspace defined by |$\mathbf {U}_A$| has as many dimensions as there are pixels in each row of the ACF array. The CCF and the ACF, however, only display a small number of modes of variability that are detectable above the noise level. Only the highest-variance components of |$\mathbf {U}_A$| are needed to capture adequately the shape changes in the ACF. The remaining low-variance components serve only to fit noise. These low-variance components may show spurious correlations with the RV signal, leading to overfitting of |$\mathbf {v}_\parallel$|⁠. It is therefore possible (and desirable) to model the activity signal adequately and avoid overfitting noise using a reduced number of dimensions.

After having repeated this procedure for all rows, we find that the reconstruction of the kth column |$\hat{\mathbf {U}}^T_k$| reproduces |$\mathbf {U}^T_{A,k}$| with good fidelity for k < 25 or so. The ratio of the median absolute deviation (MAD) of |$\mathbf {U}^T_k-\hat{\mathbf {U}}^T_k$| to |${\rm MAD}(\hat{\mathbf {U}}^T_k)$| rises to values close to unity for values of k > k_crit for which the leave-one-out cross-validation indicates that the reconstruction is poor, as shown in Fig. 6.

Following projection of the RV data into the reduced space defined by the surviving columns of |$\mathbf {\hat{U}}$|⁠, we find that the quality of the fit between the RV data and the shape model |$\mathbf {v}_\parallel$| improves rapidly at first, reaches a minimum then increases gradually as more velocity components are added to the model and overfitting starts to degrade the solution. In other words, we need even fewer principal components to model |$\mathbf {v}_\parallel$| than we need to reproduce the ACF itself.

SVD orders the principal components of the ACF in descending order of their eigenvalues |$\mathbf {S}_A$|⁠. This ordering does not take the velocity projection into account, so the ordering of principal components does not reflect accurately their contributions to the RV. Instead, we reorder the columns of |$\mathbf {\hat{U}}$| into the sequence that gives the fastest decrease in χ², obtaining the optimal fit to the RV data with the smallest number of basis vectors, as shown in Fig. 7.

Davis et al. (2017) found that 4 or 5 principal components were sufficient to capture the temporal behaviour of synthetic spectra produced by a noise-free SOAP2.0 (Dumusque et al. 2014) simulation of star-spot activity and facular suppression of convective blueshift on a rotating stellar model. As we have seen, the HARPS-N solar data contain additional profile distortions arising from changes in the instrument and Earth’s orbital motion, so more principal components are needed.

We find a good compromise between outlier clipping, number of surviving days of observation, and minimal number of basis vectors when we clip |$\mathbf {U}_A$| at 6 times the MAD, as noted above. With δχ² reordering, the χ² of |$\mathbf {v}_\perp$|is minimized at k_min = 13. We therefore use the 13 leading principal components of |$\mathbf {\hat{U}}_A$|⁠, ordered by δχ², to define the null space. We note, however, that the results that follow are only very weakly sensitive to the number of principal components used over the range 6 < k < 13 or so.

The projections of the observed velocities on to the components of |$\mathbf {U}_A$| corresponding to the 6 largest values of |$\hat{\mathbf {\alpha }}$| are plotted against time in Fig. 8.

3.4 scalpels analysis of solar data

We refer to the projection of the observed velocities on to the orthogonal complement of the time-domain scores |$\mathbf {U}_A$| of the ACF together with the outlier clipping (Section 3.2) and rank reduction (Section 3.3) algorithms collectively as Self-Correlation Analysis of Line Profiles for Extracting Low-amplitude Shifts (scalpels). The reader is referred to Appendix C1 for a concise listing of the main steps of the algorithm.

In a blind RV survey, planet-candidate detection is typically conducted using analysis of some form of periodogram such as Lomb–Scargle (e.g. Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009) or marginalized posterior versus orbital period (e.g. Mortier et al. 2015) computed from the RVs. Periodogram peaks are identified and fitted with Keplerian orbit models. This method is, however, susceptible to confusion with rotationally modulated signals from the host star.

To assess the effectiveness of the scalpels algorithm for suppressing stellar noise, we apply it to the daily averaged solar RVs in the heliocentric frame. In the absence of any dynamical shifts or instrumental calibration drift, the measured RVs should show only variations caused by line-profile shape changes arising from solar activity, changes in the instrumental point-spread function, and changes in the apparent solar rotation rate arising from Earth’s orbital eccentricity and the solar obliquity. Given that the heliocentric solar data set contains no planet signals, a frequency search for candidate periodic signals provides a means of establishing the planet-detection threshold for comparable data sets.

In Fig. 9 we show the observed heliocentric RVs minus their own mean, together with the shape-driven velocities obtained by scalpels projection and the shift-driven velocity difference between the two. The histograms of the observed velocities and the shift-driven velocities are also shown. The distribution of observed velocities is severely non-Gaussian, with a bimodal character arising from short-term (days–weeks) and long-term activity (years) variability.

After subtracting the scalpels-identified shape-driven velocity residual velocities |$\mathbf {v}_\parallel$|⁠, the shift-driven velocities |$\mathbf {v}_\perp$| are nearly constant with respect to time, with a local RMS scatter of 1.25 m s⁻¹. The Anderson–Darling 1-sample test (Scholz & Stephens 1987) indicates that the distribution of the shift-driven velocities is indistinguishable from a normal distribution (Fig. 9, lower panels), with standard deviation σ = 1.25 m s⁻¹.

Any stellar activity signature remaining in the shift-driven velocity time series is likely to show temporal correlations and departures from uncorrelated Gaussian noise. There appear to be weakly correlated residuals with amplitudes of a few tens of cm s⁻¹ on a range of time-scales upward of about 200 d. The origin of these slow drifts is unclear. They could be a shift-like manifestation of solar activity. Secular drifts induced by the instrument are also a possibility.

The Ljung–Box Q test (Ljung & Box 1978) suggests that the shift-driven velocities remain weakly correlated at all autocorrelation lags up to at least 100 d of observation. Therefore, it is likely that some activity-driven velocity components remain in |$\mathbf {v}_\perp$|⁠, but as the upper-right panel of Fig. 9 shows, they are substantially reduced relative to the original, observed velocity time series. The shape-driven signals are, as expected, strongly correlated. This offers improved detection prospects for small planetary-orbit signals at periods of tens to hundreds of d.

In Fig. 10 we show periodograms (in terms of the best-fitting semi-amplitude of a sinusoid as a function of its period) for RVs measured with the data-reduction system, transformed to the heliocentric reference frame. The periodogram of the raw velocities shows numerous candidate signals with semi-amplitudes of order 0.4 m s⁻¹, particularly between 13 and 26 d, close to the solar synodic rotation period and its first harmonic. The scalpels projection shows a very similar pattern of semi-amplitudes.

These peaks are strongly suppressed in the semi-amplitude periodogram of the shift-driven RVs (Fig. 10, bottom trace), which shows no strong frequency structure. This is important, since the background level of the periodogram peaks in the cleaned time series with no planets present, effectively sets the sensitivity for detecting planets after applying scalpels to clean the velocity measurements. Peaks with amplitudes greater than 0.30 m s⁻¹ are seen in the shift-driven RVs at P = 191.47 and 30.45 d. Their amplitudes are reduced to 0.320 and 0.339 m s⁻¹ respectively, nearly a factor of two less than those found in the observed velocities measured with the data-reduction system. We note that the excess of power at around 200 d is commensurate with the average interval between cryostat warm-ups.

4 ALGORITHM TESTS WITH PLANETS INJECTED INTO SOLAR OBSERVATIONS

We now turn to the problem of determining the impact of the scalpels signal separation on detection thresholds when weak planetary signals are present. We begin by injecting four periodic shift signals into the heliocentric CCF time series, using equation (1) to shift the rows by the small amounts required. The periods of these signals were well-spaced in log period, at non-integer periods of 7.142, 27.123, 101.543, and 213.593 d. The 27.1-d period was chosen deliberately to be close to the solar synodic rotation period. The injected signals are sinusoidal in form, with semi-amplitudes K = 40 cm s⁻¹. This is less than the amplitude of the strongest signals arising from solar variability in the upper trace of Fig. 10, but greater than the activity-driven signals remaining after subtracting the scalpels projection from the RV measurements. For a 1 M_⊙ star the corresponding planet masses are 1.2, 1.9, 2.9, and 3.7 M_⊕.

Before proceeding, we must consider the methodology used to extract velocities from the shifted CCFs and to estimate their precision from the covariances in the rows and columns of the CCF. The reader is referred to Appendix  B for the details of this methodology.

4.1 Recovery of weak injected planet signals

Following signal injection, RVs were again measured from the CCF time series using equation (B1).

Fig. 11 shows the periodograms obtained from these velocities, from the scalpels projection of the shape-driven velocity component, and from the differences between them representing pure shifts.

The periodogram of the velocities measured from the shifted CCFs does not enable us to distinguish clearly between the injected signals and RV variability intrinsic to the Sun or the instrument. The injected signals at 7.1, 27, and 100 d are detected fairly unambiguously, but the 213-d signal is suppressed and there are also many false detections of amplitude comparable to the injected signals.

The periods, semi-amplitudes and uncertainties of the five strongest signals recovered from the periodogram of the shift-driven velocities after subtraction of the shape-driven model are listed in columns 3, 4, and 5 of Table 1. Among these, four of the five strongest signals are very close to the frequencies of the injected planet signals. The mean of their semi-amplitudes is 0.428 ± 0.01 m s⁻¹, somewhat above the expected sample uncertainty of the injected values. They dominate over all residual variability and zero-point jitter signals except for a spurious 0.435 m s⁻¹ signal at P = 185.1 d. This latter period is so close to half a year that it is likely to arise from an as-yet unidentified effect of observing the Sun from the Earth, which would not be expected to affect exoplanet searches.

4.2 Simultaneous modelling of stellar variability and planetary motion

The amplitudes of the injected planet signals do not appear to be strongly attenuated in the upper trace of Fig. 11, but they are buried in a forest of activity-related peaks of similar amplitude. In the bottom trace they stand out above the suppressed activity signals. Their amplitudes could none the less be affected by activity if the planet signals themselves contaminate the scalpels projection. This could occur if the injected shift signals are not perfectly orthogonal to all elements of |$\mathbf {U}_A$|⁠, and hence partly absorbed in the scalpels projection. The data set has a finite length, so irregularly sampled superpositions of Keplerian signals will not be perfectly orthogonal to any randomly chosen vector in the same space. Moreover, a periodogram fits only a single sinusoid per frequency sample, so that cross-talk between multiple signals can lead to incorrect amplitude estimates.

The orbital perturbations of any planet and the scalpels projection process must therefore be modelled self-consistently for the signal separation to recover their semi-amplitudes as reliably as possible. Once the periods of candidate signals have been determined – either through prior knowledge of transits or via periodogram search(es) – parameter estimation and signal separation can be achieved in a single linear calculation.

The complete model of the RV data is then the sum of the model orbital velocity variations and the shape-driven velocity variations. The only unknowns are the amplitudes and phases of the orbital basis functions |$\mathbf {F}$|⁠. We can solve for these using the method of least-squares

The simultaneous solution involves computing the shape-driven variations from the difference between the observed RVs and the model velocities. In the projection-operator language of Section 3, we solve for the vector |$\mathbf {\theta }_{\rm orb}$| that minimizes |$\chi ^2=(\mathbf {P}_\perp \cdot \delta \mathbf {v}^T) \cdot \Sigma ^{-1} \cdot (\mathbf {P}_\perp \cdot \delta \mathbf {v})$|⁠, where |$\delta \mathbf {v}\equiv \mathbf {v}_{\rm obs}-\left\langle v_{\rm obs}\right\rangle -\mathbf {F}\cdot \mathbf {\theta }_{\rm orb}$|⁠.

Here we use the simplifying assumption that the RV measurements are uncorrelated. The covariance matrix is diagonal and its log determinant is |$\sum _{j=1}^{l}\ln \mathbf {\Sigma }_{jj}$|⁠. The diagonal elements |$\Sigma _{jj}={\rm Var}(\mathbf {v}(t_j))$| are calculated using equations (B2) and (B8). If the RV data are sufficiently densely sampled, a time-dependent covariance model with a kernel incorporating the stellar rotation period and active-region lifetime could also be included – see, e.g. Gilbertson et al. (2020). For convenience, we summarize the algorithm in Appendix C2.

In Table 1, columns 6–8 list the amplitudes and uncertainties of the sinusoidal signals recovered from the data at the known periods of the injected signals. The standard deviation of the four individual recovered semi-amplitudes is σ = 0.010 m s⁻¹. The sample mean and standard deviation (⁠|$\sigma /\sqrt{4}$|⁠) of the signal amplitudes are 0.433 ± 0.005 m s⁻¹, again somewhat above the injected value. The four individual signals deviate from the injected values by amounts that are consistent with their individual estimated 0.066 m s⁻¹ semi-amplitude uncertainties.

The improvement in signal separation is also apparent from Fig. 12. The top two traces are almost the same as those in Fig. 11, but the balance of the signal separation is changed by the explicit modelling of the orbital motion at the known periods. The periodogram of the fitted orbital model illustrates the apparent signal attenuation that can occur when fitting multiple signals with a single sinusoid. The final residuals are very similar to those of Fig. 10.

The precision of the recovered semi-amplitudes is poorer when the velocities are left uncorrected for profile-shape variations, by setting the dimension k_max of the null space to zero. If sinusoids are fitted to the raw |$\mathbf {v}_{\rm obs}$| at the same four periods, we obtain semi-amplitudes 0.466, 0.560, 0.481, and 0.129 ms ⁻¹, whose sample mean and standard deviation are 0.409 ± 0.083 m s⁻¹. Fig. 13 shows clearly the improvement in fidelity of the recovered amplitudes when the optimal shape model is applied.

The formally propagated 0.066 m s⁻¹ uncertainties in the semi-amplitudes are |$\sim 1.6\sigma _{\rm v} /\sqrt{N_{\rm obs}}$| for N_obs = 853 if we adopt a single-measurement precision of σ_v ≃ 1.25 m s⁻¹ based on the RMS scatter in the heliocentric solar velocities after removal of the shape perturbations (see Section 3.4).The RV amplitude precision appears to scale as expected for uncorrelated random variables to a level ∼20 times better than the single-measurement precision.

The average of the recovered semi-amplitudes shows no evidence of bias relative to the injected value when signal separation is performed with prior knowledge of the orbital period, as is the case with transiting planets. Thus, the RV semi-amplitudes inferred from the cleaned velocities can be significantly more reliable than RV semi-amplitudes inferred from original velocity measurements.

We conclude that the scalpels method succeeds in reducing correlations between apparent velocities due to stellar variability, based solely on line shape changes and without making use of time-domain information. This decoupling from time-domain information allows planet signals to be recovered with good fidelity even when they fall close to the stellar rotation period, as is the case with the 27-day signal injected here.

5 DISCUSSION

5.1 Summary

We have presented a new algorithm for extracting precise RV estimates from high-resolution spectroscopic planet surveys. The algorithm begins with a list of CCFs for each observation epoch, constructs a reduced-rank representation of stellar variability and reconstructs CCFs which have been cleaned of most stellar variability. We demonstrated the algorithm using observations of the solar spectrum from HARPS-N. We verified and validated that the algorithm can accurately detect multiple simulated planets injected into solar observations, spanning a wide range of orbital periods.

5.2 Planets injected into solar observations

Our algorithm recovered the RV signals of injected planets with semi-amplitudes of just 40 cm s⁻¹ with an SNR of 6 based on a time series of 853 CCF measurements.

The semi-amplitude uncertainty of 6.6 cm s⁻¹ implies that it is possible for intensive Doppler spectroscopy campaigns to measure the masses of transiting planets (i.e. with well-measured orbital period and inclination) at |$\sim 15{{\ \rm per\ cent}}$| precision for RV semi-amplitudes as small as K = 40 cm s⁻¹ for a solar twin, even if there are multiple planets spanning a wide range of orbital periods.

The precision with which we can measure velocity semi-amplitudes is insensitive to the velocity semi-amplitude of the planet (for small velocity semi-amplitudes). Thus, our results suggest that intensive Doppler spectroscopy campaigns could detect and measure the masses (times sine of inclination angle) to 6σ precision for planets with RV semi-amplitudes of ∼60 cm s⁻¹, if there were independent evidence for a planet at a given period (e.g. transit, direct imaging) orbiting a sufficiently bright (V ≃ 6) Sun-like star.

5.3 Areas for Further Research

The algorithm presented successfully mitigates the impact of solar variability. Nevertheless, there are multiple lines for additional research that are likely to lead to further improvements in the Doppler sensitivity.

5.3.1 Simulations using realistic observing cadence

For planet injection tests in Section 4, we used CCFs from HARPS-N solar observations taken on 853 of a total 886 d. The observations span 1681 d, and over 70 per cent of the observations are spaced by just 1 d. While there are seasonal shutdowns each year and gaps of up to two weeks, this observing cadence is more favourable than that of most targets of Doppler planet surveys. We recommend that future simulations of Doppler planet surveys explore how the planet detection sensitivity depends on observing cadence in the presence of intrinsic stellar variability and after incorporating advanced methods for mitigating stellar variability such as presented here. Hall et al. (2018) have carried out such a study for the HARPS-3 project, which will produce stellar data sets of comparable duration and quality to the solar data set examined here, albeit with seasonal gaps. A similar sensitivity study using HARPS-N solar data has been published recently by Langellier et al. (2020), who concluded that a decade or more of observations could be needed to achieve a 5σ detection of a 50 cm s⁻¹ signal with a period of 225 d, given an instrumental white-noise uncertainty of 80 cm s⁻¹ and a perfect model of activity-driven RV variability. Haywood et al. (2020) analysed an 8-yr sequence of synthetic solar RVs derived from SDO images, with a six-month duty cycle and decorrelation against hemispherically averaged magnetic field. They succeeded in recovering an injected signal with K = 30 cm s⁻¹ and P = 300 d. Our results, using real data but an idealized observing pattern, give similar cause for optimism.

5.3.2 Granulation

The CCF time series used for verifying our algorithm was based on 15-min daily averages of the HARPS-N solar CCFs. Each exposure time is 300 s, yielding an SNR in the range 250 < SNR < 400. This is comparable to that of exoplanet observations for a host star of magnitude 5.5 using the same exposure time. Combining three contiguous such exposures in a 15-min block is reasonably effective at averaging-out the side effects of p-mode pulsations which occur on a ∼5 min time-scale (Chaplin et al. 2019). Modern exoplanet surveys typically choose an exposure time of at least 15 min, averaging subexposures if necessary to average out spectral variability due to pulsations and avoid saturation. Such exposures are repeated at intervals of 2–3 h to mitigate the effects of granulation. Our use of a single 15-min block per day is likely to be less effective at eliminating intrinsic stellar variability due to the granulation pattern, as studied in detail by Meunier et al. (2015). We tested this by averaging all data satisfying the data-selection constraints in Section 2 throughout each day of observations, rather than down-selecting to a single 15-min block. The longer daily baseline reduced the day-to-day scatter in |$\mathbf {V}_\perp$| from 1.25 to 1.08 m s⁻¹. This is less than the improvement expected if photon counts were the limiting noise source, but consistent with the improvement expected through averaging over velocity fluctuations arising from photospheric granulation noise.

In principle, granulation might imprint on the stellar spectrum differently than variability on the time-scales of magnetic activity (Cegla 2019). It is likely that some of the 13 modes of RV variability found in our data may be attributable to granulation.

5.3.3 Searching for low-mass planets with outer giants

Some complexities of planetary-system architecture are beyond the scope of this initial study, but merit further investigation. For example, signals from terrestrial planets with sub-m s⁻¹ amplitudes may be superposed on much larger signals of giant planetary companions at longer orbital periods. If the giant’s orbit is well-characterized, it can be included directly in the model. Indeed the barycentric data considered here contain such a signal, with the synodic period of Jupiter. Rather than simply work in the heliocentric frame, we repeated the analysis of Section 4.2 by injecting the same four signals into the barycentric CCFs. To mimic the situation where the giant’s orbit is incomplete or poorly determined, we performed a GP regression with a squared-exponential covariance kernel to smooth out variations on time-scales longer than 300 d. All four of the injected signals were recovered successfully at the correct amplitudes.

More complex cases might involve transiting systems in which a temperate Earth-sized object is accompanied by a few strongly interacting compact short-period companions. In this case, gravitational interactions might make the fully linear algebra approach of the method less effective. As shown in Fig. 9, however, signal separation and detection is reasonably effective even in a ‘blind-search’ scenario where simultaneous linear fitting is not possible. Furthermore, even a nonlinear parametric model of the orbital reflex motion can be included as a fitting function and its parameters optimized together with the coefficients of ACF basis vectors in an MCMC scheme, at the cost of some additional computational overheads.

5.3.4 Reconstructing full spectra

In this paper, we reconstructed the shape-driven velocity signal using a reduced-rank representation and scores derived from the ACF of the CCFs. In principle, the same approach could be used to reconstruct the CCFs themselves or even the full spectra. Then, the reconstructed spectra could become objects for further analysis, e.g. performing line-by-line analysis (Dumusque 2018). We anticipate that a reconstruction of the spectra would likely lead to more significant deviations from the observed spectra than those we found when reconstructing the velocities alone from the ACF. This is because the CCF, from which the ACF is derived, combines information from a large list of spectral lines, even though different lines are known to respond differently to stellar variability (Dravins et al. 1981; Toner & Gray 1988). Furthermore, variability of the continuum contributes very little to the CCF. Future studies could analyse the residual spectra to gain insights into more subtle ways in which stellar variability manifests in the observed spectra (e.g. Davis et al. 2017; Thompson et al. 2017; Dumusque 2018; Wise et al. 2018; Cretignier et al. 2020).

5.3.5 Multiple CCF masks

The results of full-spectrum reconstruction could inform the development of additional stellar variability indicators and/or the design of CCF masks. As noted above, summarizing the spectrum with a single CCF averages over the different responses of lines which are more or less sensitive to stellar activity. For this study, we have used the CCF generated by the HARPS-N DRS. However, one could design multiple CCF masks to compute multiple CCFs at each observation epoch. The choice of mask could be based on astrophysical insights or machine-learning approaches. Either way, each mask would include spectral lines exhibiting common patterns of response to stellar variability, so as to provide greater sensitivity for recognizing line shape variations.

Once CCFs have been computed using multiple masks, the ACF can be computed for each CCF-mask pair. The underlying algorithm presented in Section 3 can be applied to compute time-domain subspaces either separately for each mask or simultaneously. The resulting velocities derived from each mask could be analysed separately for diagnostic purposes and then combined for inference (e.g. Zechmeister et al. 2018). The use of multiple masks naturally leads to multivariate time series for estimates of the velocity and stellar variability indices.

5.3.6 Time-domain modelling

Stellar variability generally distinguishes itself from planet-induced orbital motion both in the wavelength domain (e.g. relative depths of lines, line shapes) and in the time domain (e.g. deviations from a strictly Keplerian signal). We caution, however, that subtler effects such as activity-driven changes in the stellar gravitational redshift may manifest as pure shifts at amplitudes up to 10 cm s⁻¹ (Cegla et al. 2012). In this paper, we have described an improved algorithm for measuring RVs in the presence of stellar variability that does not make use of time-domain information. That is, the computation of |$\mathbf {v}_\parallel$| depends only on the ensemble of measured CCFs, but not on the times at which the observations were taken. For the purpose of demonstrating our algorithm, we have used traditional maximum likelihood estimation based solely on the cleaned velocities. Our algorithm naturally produces additional stellar variability indicators [i.e. |$\mathbf {U}_A(t)$| and |$\mathbf {U}_C(t)$|] that could be modelled simultaneously with the velocities, following methods developed by Rajpaul et al. (2015) and generalized in Jones et al. (2017).

We recommend that future simulations explore the potential for further improvements in the sensitivity to small planets by combining algorithms for utilizing information in the wavelength and time domains.

5.3.7 Integrating into Doppler planet survey toolboxes

The algorithm developed in his paper can be implemented efficiently using a standard linear algebra toolbox. Thus, it can be readily integrated into existing or future software packages to analyse Doppler planet search observations. For example, one could inspect the posterior marginalized over all parameters except orbital period (i.e. a periodogram), similar to Mortier et al. (2015), but replacing the standard likelihood with equation (14). When considering observations of a star potentially hosting multiple planets, a periodogram of the reconstructed CCFs could be applied iteratively, removing one signal at a time. Alternatively, one could apply sparse regression techniques, i.e. fit for the semi-amplitudes of many periodic signals simultaneously, while applying regularization to the semi-amplitude, so as to find a family of maximum likelihood solutions for each plausible number of planets.

This could be implemented efficiently for an a priori unknown number of signals using either the alternating direction method of multipliers (ADMM) algorithm or the spectral projected-gradient algorithm (as in Hara et al. 2017). Once the putative orbital periods have been identified, then MCMC-based techniques can be used to perform parameter estimation (e.g. Ford 2006; Nelson, Ford & Payne 2014) replacing the standard likelihood with our equation(14). Then, one could compare the Bayes factor or ratio of marginalized likelihoods (i.e. ‘evidences’) for models with various numbers of planets. While there is still active research in finding efficient and robust algorithms for estimating Bayes factors for models with several planets, replacing the likelihood with equation (14) would be algorithmically straightforward and is expected to add only a modest additional cost.

5.3.8 Design of Doppler planet surveys

The observing strategy of Doppler planet surveys (e.g. number of observations per star, distribution of duration between observations) can have a significant impact on the sensitivity for detecting planets and its dependence on orbital period. Simulations of Doppler planet surveys have been used to inform survey design choices (e.g. Ford 2008; Burt et al. 2018; Cloutier et al. 2018; Hall et al. 2018).

Previous studies have typically ignored or adopted simplistic models for intrinsic stellar variability (e.g. assuming that a fixed fraction of spurious velocities due to stellar variability can be corrected). Our algorithm for simultaneously inferring the effects of stellar variability and planetary signals could be readily incorporated into more advanced survey simulations incorporating explicit models of the effects of stellar activity such as that published recently by Damasso et al. (2019).

Now that there is a concrete and computationally efficient strategy for mitigating stellar variability, we recommend that future studies conduct new survey simulations to compare candidate Doppler planet survey strategies. We anticipate that our method will significantly increase the value of observing the same star many times, since we substantially reduce the correlation of inferred velocities at different times. On the other hand, it may be that the improved ability to mitigate stellar activity reduces the importance of obtaining a dense set of observations over the rotation period or active region lifetime. Monte Carlo simulations are necessary to understand the interaction of these two effects and the implications for future planet surveys. In addition to informing survey strategies that do not make use of information from previous observations, our algorithms could be folded into adaptive scheduling algorithms that maximize a merit function (or minimize a cost function), so as to maximize the efficiency of a Doppler survey (e.g. Ford 2008).

ACKNOWLEDGEMENTS

This article is a result of collaborative scholarly efforts from the residency of the Research Group on Big Data and Planets at the Israel Institute for Advanced Studies. We acknowledge valuable discussions with Tsevi Mazeh, Eric Feigelson, Shay Zucker, Lev Tal-Or, and Dovi Poznanski. We thank the anonymous referee for numerous insightful suggestions that led to major improvements, including the use of leave-one-out cross-validation and streamlining of the mathematical approach used in Section 4.2. ACC acknowledges support from the Science and Technology Facilities Council (STFC) consolidated grant number ST/R000824/1 and UKSA grant ST/R003203/1. EBF acknowledges support from NSF award #1616086, NASA grant #80NSSC18K0443, and Heising-Simons Foundation grant #2018-0851. This work was supported by a grant from the Simons Foundation/SFARI (675601, EBF). EBF acknowledges the support of the Ambrose Monell Foundation and the Institute for Advanced Study. EBF acknowledges support from the Penn State Eberly College of Science, Department of Astronomy & Astrophysics, Institute for Computational & Data Sciences, the Center for Exoplanets and Habitable Worlds, the Center for Astrostatistics. DFP acknowledges NASA award number NNX16AD42G. This work was performed partly under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute (RDH). AM acknowledges support from the senior Kavli Institute Fellowships, funded by the Kavli Foundation. XD is grateful to the Branco Weiss Fellowship – Society in Science for its financial support. HMC acknowledges financial support from the National Centre for Competence in Research (NCCR) PlanetS, supported by the Swiss National Science Foundation (SNSF), as well as a UK Research and Innovation Future Leaders Fellowship. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Framework Programme (grant agreement No. 851555). This work has been carried out in the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). The authors gratefully acknowledge the support of the International Team 453 by the International Space Science Institute (Bern, Switzerland). Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain). The HARPS-N project has been funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origins of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), the Italian National Astrophysical Institute (INAF), the University of St Andrews, Queen’s University Belfast, and the University of Edinburgh. The citations in this paper have made use of NASA’s Astrophysics Data System Bibliographic Services. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch.

DATA AVAILABILITY

The HARPS-N solar data products for the first 3 yr of the period described in this paper are publicly available through the DACE platform (https://dace.unige.ch) developed in the frame of PlanetS. All other research data underpinning this publication (https://doi.org/10.17630/9ec51274-cbea-445f-b188-112f4734c6e9) and the python code and notebook used to prepare all diagrams in this paper (https://doi.org/10.17630/957d6320-1969-43c8-b6bb-a8ae53d1f657) will be made available through the University of St Andrews Research Portal.

REFERENCES

APPENDIX A: CALCULATION OF PROFILE DERIVATIVES

Derivatives of the rows of the CCF with respect to velocity are required for applying RV shifts to the rows of the CCF in both correction from the barycentric to the heliocentric reference frame (Section 2.1) and for injecting synthetic orbital reflex motion signals (Section 4).

Simple numerical derivatives computed directly from the data have the undesirable property that they amplify noise. If the profile shape was time-invariant, this problem could be overcome by fitting the Taylor-series model derived from the mean profile, to a sequence of CCFs shifted by orbital reflex motion. However, magnetically active regions rotating across the star’s visible hemisphere distort the CCF profile (Meunier et al. 2010). This time variation in the shape of the profile alters the RV measured by fitting a fixed profile or parametric function to the CCF. Currently, it is unclear whether the cross-terms between such intrinsic stellar variability and small Doppler shifts caused by planets are sufficiently small to be neglected. Therefore, this paper attempts to estimate the derivative of the CCF at each observation time and focuses on solar data for which very high signal-to-noise data is available.

Rather than estimating the derivative from a constant mean spectrum, we estimate the derivative of the CCF at each epoch t_j based on a reduced-rank reconstruction of the CCF using the mean profile and the k_max = 10 leading principal components obtained by singular-value decomposition (SVD) of the full ensemble of CCFs. The optimal choice of k_max that best separates the effects of profile shape changes from the effects of oversampled photon noise is calculated as described in Appendix B1.

APPENDIX B: TAYLOR-SERIES VELOCITY MEASUREMENT FROM CCF

In order to ensure that the sampling of the CCF did not bias the scale of the derived velocities, we performed linear regression of both the DRS velocities and the Taylor-series velocities against the solar barycentric RV computed with JPL horizons (Giorgini et al. 1996). The scale factors differed from unity by 0.8 and 2.4 per cent, respectively. The barycentric RVs were scaled by these factors for all subsequent calculations, to ensure that the velocities were transformed correctly to the heliocentric reference frame. The agreement between the two sets of corrected velocities is illustrated in Fig. B2.

B1 Covariance matrix for time-varying CCFs

For calculating instantaneous RVs and their precisions from equations (B1) and (B2), the covariance matrix Σ must contain information about the SNR of each observation and the correlations between neighbouring pixels. The covariances between different pixels in the time series of CCFs fall into two categories.

The first is systematic uncertainty arising from temporal variability of the profile shape. The sample covariance matrix |$\Sigma = {\rm Cov}(\mathbf {R}) = \mathbf {R}^T\cdot \mathbf {R}/m$|estimates the systematic covariances between the columns of |$\mathbf {R}$|⁠, where each column has m elements, one per row of the CCF time series. These arise from correlations between different parts of the line profile.

The second is random measurement error arising from the finite SNR of the original spectrum. The CCF sampling interval is matched to the velocity increment per CCD pixel of the instrument. Although neighbouring pixel values in the original spectrum are statistically independent, interpolation during rebinning and calculation of the CCF introduces correlations between neighbouring CCF samples. We therefore need to account for the correlations between CCF elements at neighbouring values of velocity.

The dispersion of HARPS-N gives a near-constant velocity increment per physical CCD pixel in the dispersion direction of 0.82 km s⁻¹ (Dumusque et al. 2020). The DRS delivers a default CCF that is also sampled in velocity increments of 0.82 km s⁻¹, so uncorrelated photon-noise fluctuations in the extracted spectra are spread over about 2 CCF velocity increments after rebinning and interpolation.

The systematic covariances and the instantaneous correlated-noise pattern are unrelated to each other, and enter into the velocity and error calculation in different ways. To estimate the independent systematic uncertainty that affects every observation through profile shape changes, the systematic covariance matrix should be used with equation (B2), separately from the small-scale noise pattern.

When this reduced-rank version is subtracted from the full covariance matrix, we obtain the covariance matrix of the remaining high-frequency noise, as shown in the middle panel of Fig. B1. Its strongest feature is a ridge of covariance with an approximately triangular cross-section, whose amplitude depends on the average SNR of the original spectra, and whose width reflects the sampling of the CCF. The peak variance along its diagonal is about 100 times smaller than the peak amplitude of the large-scale covariance pattern.

We optimize k_max to give a clear separation between the large-scale column covariances and the covariances between neighbouring elements in each row arising from oversampling of photon noise. To achieve this we divide the peak value in the bottom panel of Fig. B1 by the range of all other diagonal means with lags differing by more than 1 km s⁻¹ (slightly more than the CCF sampling interval) in velocity from the leading diagonal. This ratio shows a well-defined peak at k_max = 6, where the separation between the two variance patterns is optimized.

B2 Comparison with DRS velocities

We use this covariance matrix with equations (B1) and (B2) to estimate the RVs |$\mathbf {\delta v}(t_j)$| and their variances due to photon noise and profile-shape changes.

In Fig. B2 we show that RVs measured using this method are almost identical to those reported by the HARPS-N DRS. The overall RMS scatter in the differences between the two estimates of the velocity is 0.044 m s⁻¹.

The formal errors derived from equation (B2) are typically 0.11 m s⁻¹, though the data set includes a small number of days of sparse data yielding uncertainties greater than 0.2 m s⁻¹. Such small error estimates are not surprising. The formal uncertainties computed for individual observations by the DRS indicate an average of 0.24 m s⁻¹ photon noise per exposure. We average many such exposures per day, so the photon noise is insignificant in comparison to the systematic uncertainties arising from profile-shape variations, zero-point errors and calibration drift.

The independent systematic variances arising from large-scale profile-shape changes are re-computed from equation (B2) using the reduced-rank covariance matrix |$\Sigma = {\rm Cov}(\mathbf {R}_{k_{\rm max}})$| from equations (B3) and (B4). The resulting systematic error per observation is 0.82 ± 0.017 m s⁻¹, which is very close to the RMS amplitude of the scatter in Fig. 9.

APPENDIX C: scalpels ALGORITHM

C1 Simple blind search

Input:|$\mathbf {C}$| (2D CCF time series), |$\mathbf {x}$| (CCF array velocities), |$\mathbf {v}_{\rm obs}$| (RV measurements), |$\sigma (\mathbf {v}_{\rm obs})$| (RV error estimates).

• Compute ACF with equation (3) and normalize.

• Compute SV decomposition of ACF with equation (5).

• Perform MAD clip at appropriate threshold on columns of |$\mathbf {U}_A$|⁠.

• Recompute SVD of surviving ACF rows.

• Withhold each row of ACF in turn, perform SVD, reconstruct corresponding row of |$\mathbf {\hat{U}}_A$| (equations 9, 10).

• Keep columns of |$\mathbf {U}_A$| s.t. |${\rm MAD}(\mathbf {U}^T_k-\hat{\mathbf {U}}^T_k)/{\rm MAD}(\hat{\mathbf {U}}^T_k)\ll 1$|⁠.

• Compute |$\hat{\mathbf {\alpha }}$| (equation 6).

• Rearrange elements of |$\hat{\mathbf {\alpha }}$| and columns of |$\mathbf {U}_A$| in order of descending |$|\delta \chi ^2(\mathbf {v}_\perp)|$| (Sect 3.3).

• Select number of basis vectors to keep (Sect 3.3).

• Compute |$\mathbf {v}_\parallel$| (equation 7) in reduced-rank basis.

• Compute |$\mathbf {v}_\perp =\mathbf {v}_{\rm obs}-\mathbf {v}_\parallel$|⁠.

Return:|$\mathbf {v}_\parallel$| (shape-driven RV), |$\mathbf {v}_\perp$| (shift-driven RV).

C2 Simultaneous sinusoidal fit

Input:|$\mathbf {C}$| (2D CCF time series), |$\mathbf {x}$| (CCF array velocities), |$\mathbf {v}_{\rm obs}$| (RV measurements), |$\sigma (\mathbf {v}_{\rm obs})$| (RV error estimates), {ω₁, ⋅⋅⋅, ω_n} (orbital frequencies for n planets).

Perform steps (i) - (vii) above and see Section 4.2.

• Compute |$\mathbf {F} = \lbrace \cos \omega _1 t_j, \sin \omega _1 t_j, \cdots , \cos \omega _n t_j, \sin \omega _n t_j\rbrace$|⁠.

• Compute reduced-rank covariance matrix (equation B5).

• Compute row variances |$\sigma ^2_j$| of |$\mathbf {R}-\mathbf {R}_{k_{\rm max}}$| (Appendix  B).

• Compute model of full covariance matrix (equation B8).

• Compute |$\mathbf {C^{\prime }}$| and |${\rm Var}(\mathbf {v}(t_j))$| (equation B2).

• Construct |$\Sigma = {\rm Diag}({\rm Var}(\mathbf {v}(t_j)))$|⁠.

• Construct |$\mathbf {P}_\perp =(\mathbf {I}-\mathbf {U}_A\cdot \mathbf {U}_A^T)$| in reduced-rank basis.

• Compute |$\mathbf {v}_\perp = \mathbf {P}_\perp \cdot \mathbf {v}_{\rm obs}$|⁠.

• Compute |$\mathbf {F}_\perp = \mathbf {P}_\perp \cdot \mathbf {F}$|⁠.

• Solve equation (13) to obtain |$\mathbf {\theta }_{\rm orb}$|⁠.

• Compute |${\rm Var}(\mathbf {\theta }_{\rm orb})=1/{\rm Diag}(\mathbf {F}_\perp ^T \cdot \mathbf {\Sigma }^{-1}\cdot \mathbf {F}_\perp)$|⁠.

• Compute |$\mathbf {v}_{\rm orb} = \mathbf {F}\cdot \mathbf {\theta }_{\rm orb}$|⁠.

• Compute |$\mathbf {v}_{\rm resid} = \mathbf {v}_\perp -\mathbf {F}_\perp .\theta _{\rm orb})$|⁠.

• Compute |$\mathbf {v}_\parallel =\mathbf {v}_{\rm obs}-\mathbf {v}_\perp$|⁠.

Return: RV amplitudes and variances, |$\mathbf {v}_\parallel$|⁠, |$\mathbf {v}_\perp$|⁠.

Author notes