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

Stellar magnetic activity produces time-varying distortions in the photospheric line proﬁles 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 identiﬁers and spanned by the amplitude coefﬁcients 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 − 1 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 − 1 , 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.


I N T RO D U C T I O N
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 highresolutionéchelle spectrographs have allowed exoplanet surveys to probe planets with radial velocity (RV) semi-amplitudes of just ∼1 m s −1 . 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 E-mail: acc4@st-andrews.ac.uk † NASA Sagan Fellow. 30 cm s −1 precision per RV observation on Proxima Cen (Suárez Mascareño et al. 2020), and EXPRESS has achieved 58 cm s −1 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.

C RO S S -C O R R E L AT I O N F U N C T I O N A N D I T S AU TO C O R R E L AT I O N
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 −1 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 ) 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 principalcomponent 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 ) = C(v i , t j ) − C(v i ) obtained by subtracting the time-averaged profile of the entire 5yr 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 −1 . In summer these conditions are satisfied for up to 4 h d −1 , and in winter for up to 2 h d −1 . Days on which no data were obtained are not shown, so although time increases up the vertical axis, the time-scale is not linear. The barycentric residual CCF (upper left) shows the solar reflex motion about the Solar system barycentre. Dominated by Jupiter, it has a semi-amplitude of 12 m s −1 at Jupiter's synodic period of 398 d. The residual CCF in the heliocentric frame (lower left) shows a secular change in line depth, deepening as the solar activity level declines over the period of observation. The residual ACF (right panels) shows temporal variability that is correlated with that of the CCF (compare left and right panels), but is unchanged by Doppler shifts being applied to the CCF (compare upper right and lower right panels).
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 signalto-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.

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).
We use a Taylor-series approximation to eliminate the solar barycentric motion from the CCF time series, adding scaled derivatives of the instantaneous profile shape at time t j to the barycentric CCF: 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 −1 , which is much less than the h = 820 m s −1 sampling interval between neighbouring CCF elements. The truncation error in equation (1) is therefore significantly less than 3 /12h 3 2.6 × 10 −8 .
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 −1 . 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 −1 .
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.

Autocorrelation of CCF and shift invariance
Since our goal is to separate the effects of genuine dynamical Doppler shifts from spurious shifts caused by line-shape changes, we aim to characterize changes in CCF profile shape in a way that is invariant to translation in velocity space. The autocorrelation function (ACF) of the CCF has the desired property that it is invariant to translation (Adler & Konheim 1962). The ACF A(δv) is the expectation value of the vector cross-product of the CCF with itself at a sequence of lags δv: For the sake of brevity, in this manuscript, we refer to the ACF of the CCF as simply 'the ACF'.
The CCF series has m rows representing individual observations and l columns representing individual velocity bins. We compute the ACF of every CCF in the time series. This is done by sequentially shifting the CCF by integer numbers of velocity steps, or CCF 'pixels', modulo the number l of elements in the CCF, and comultiplying by the unshifted CCF: This set of circular shifts and cross-products is repeated for every observation, to obtain a time sequence of ACFs. The m rows of the ACF time series have the same length l as the original CCFs and are normalized to a mean value of unity. There is sufficient pseudo-continuum to either side of the dip in the CCFs to ensure that this circular autocorrelation procedure is sensitive to long-range correlations while minimising edge effects. 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.

Principal-component analysis of stellar variability
The principal modes of variability in the CCF can be isolated by calculating the singular-value decomposition (SVD) of the ensemble of CCFs: The same method yields the principal components of the ensemble of ACFs: The diagonal matrices S C and S A list the singular values (eigenvalues) of the principal components in decreasing order. Fig. 3 (top) shows the eigenvectors P C,k (v i ) and P A,k (v i ) (also known as loadings) of the leading (k = 1···6) principal components of the heliocentric CCF (left) and ACF (right) time series. They represent orthonormal modes of profile variability. The columns of U C,k (t j ) and 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 continuousflow 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 profileshape 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.

T H E RV R E S P O N S E TO AC F T I M E -D O M A I N VA R I AT I O N S
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.

Projection into the ACF time-domain subspace
The set of RV observations has m elements, and can be thought of as a vector v obs belonging to an m-dimensional space S. The l We first subtract the inverse-variance weighted mean v obs from the vector v obs of RVs measured with the data-reduction pipeline, to ensure orthogonality. We project the difference v obs − v obs on to the time-domain manifold spanned by the basis formed by the columns of the matrix U A . Each row of U A corresponds to an individual observation identified with a unique time-stamp. For conciseness we refer to this observation-identifier space as the 'time domain'. The inner product of the kth column of U A (i.e. the scores associated with kth basis vector for ACF decomposition) with the velocities yields the response α k = U T A,k · (v obs − v obs ) of the RV to the time variation of U A,k . The vector of response factors is then The sum of the scaled velocity contributions from all principal components of the ACF is then v = U A ·α. This velocity vector v lies within the subspace U, and gives a complete model of the RV perturbations arising from the changes in profile shape to which the ACF is sensitive: The product P = U A · U T A is a projection operator. The operator lie outside the subspace U A , and therefore v ⊥ preserves the shifts to which the ACF is insensitive. Information is, however, lost in this process. The resulting velocities v ⊥ are biased down by the projection of the velocities on to U A , as discussed in Section 4. From here on, we refer to v obs as the 'measured' or 'observed' velocities'; v as the 'model' or 'shape-driven' velocities; and v ⊥ as the 'shift-driven' velocities. Fig. 5 shows that the shapedriven velocity perturbations, v , are strongly correlated with the observed velocities, reproducing faithfully the long-term and shortterm fluctuations dominated by stellar activity.
Given a set of measured RVs (v 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 shapedriven perturbations to the RV (v ), so as to provide a substantially cleaned set of shift-driven RVs (v ⊥ ).

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 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 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.

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 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 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 v . It is therefore possible (and desirable) to model the activity signal adequately and avoid overfitting noise using a reduced number of dimensions.
To determine the optimal size of the null space, we used leaveone-out cross-validation (Celisse 2014). Holding out each row A j of the ACF in turn, we decompose the remaining rows and compute the singular-value decomposition: We reconstruct an estimateÛ j of the missing jth row of U A by fitting the eigenvectors and eigenvalues to the jth row of the ACF: After having repeated this procedure for all rows, we find that the reconstruction of the kth columnÛ T k reproduces U T A,k with good fidelity for k < 25 or so. The ratio of the median absolute deviation (MAD) of U T k −Û T k to MAD(Û 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.  . As more dimensions are added to the subspace defined by the time-domain coefficientsÛ of the principal components of the ACF, the χ 2 of v ⊥ decreases to a minimum value, then increases gradually as overfitting degrades the solution. The minimum is reached more rapidly when the columns of U are sorted in order of decreasing δχ 2 (blue) rather than in order of their corresponding singular values (orange). With δχ 2 sorting, the optimal size of the null space is defined by the minimum at k = 13.
Following projection of the RV data into the reduced space defined by the surviving columns ofÛ, we find that the quality of the fit between the RV data and the shape model v 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 v than we need to reproduce the ACF itself.
SVD orders the principal components of the ACF in descending order of their eigenvalues 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Û into the sequence that gives the fastest decrease in χ 2 , 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 U A at 6 times the MAD, as noted above. With δχ 2 reordering, the χ 2 of v ⊥ is minimized at k min = 13. We therefore use the 13 leading principal components ofÛ A , ordered by δχ 2 , 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 U A corresponding to the 6 largest values ofα are plotted against time in Fig. 8.

SCALPELS analysis of solar data
We refer to the projection of the observed velocities on to the orthogonal complement of the time-domain scores 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 shiftdriven 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 v , the shift-driven velocities v ⊥ are nearly constant with respect to time, with a local RMS scatter of 1.25 m s −1 . 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 −1 .
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 −1 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 The Ljung-Box Q test (Ljung & Box 1978) suggests that the shiftdriven 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 v ⊥ , but as the upperright 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 semiamplitude 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 −1 , 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 −1 are seen in the shiftdriven RVs at P = 191.47 and 30.45 d. Their amplitudes are reduced to 0.320 and 0.339 m s −1 respectively, nearly a factor of two less than those found in the observed velocities measured with the datareduction system. We note that the excess of power at around 200 d is commensurate with the average interval between cryostat warm-ups.

A L G O R I T H M T E S T S W I T H P L A N E T S I N J E C T E D I N TO S O L A R O B S E RVAT I O N S
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 −1 . 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. Figure 9. Upper left: observed RVs transformed to the heliocentric reference frame, together with their SCALPELS-separated shape-driven and shift-driven velocity components, offset by ±10 m s −1 for clarity. Upper right: ACFs of the RVs, corrected for missing dates of observation and normalized to unity at zero lag. The original RVs (blue) are seen to be strongly correlated at all time lags, and the shape-driven velocities (orange) even more so. The autocorrelation of the shift-driven velocities decays substantially more rapidly than the observed or shape-driven velocities. Any correlation at lags longer than ∼70 d is negligible. The red curve is the ACF of the shift-driven velocities shuffled into random order, and is effectively uncorrelated. Lower left: Histograms show that the RMS scatter has been reduced from 1.76 m s −1 in the original data to 1.25 m s −1 in the shift-driven velocities. Lower right: shift-driven velocities sorted in ascending order and overplotted with the cumulative normal distribution with σ = 1.25 m s −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 −1 , 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 −1 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.

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 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 Figure 10. Semi-amplitude periodogram demonstrating the value of signal separation by projection of the observed RVs on to the principal timedomain components of the ACF. The top periodogram (blue) is for the measured velocities derived from the CCF. The second trace (orange) is for the shape-driven velocities v produced by the SCALPELS projection. The third periodogram (green) is that of the shift-driven velocities v ⊥ remaining after subtraction of the shape-driven velocities from the observations. Lightblue bars denote the approximate ranges of the solar rotation period and its first harmonic, and periods of 6 months and 1 yr. Figure 11. Periodograms of velocities derived from the heliocentric solar CCFs when four sinusoidal signals of 40 cm s −1 have been injected at the periods denoted by the vertical blue lines. The traces are as defined in the caption of Fig. 10. The uncertainty in the amplitude of the sinusoid is almost independent of period. Its 1σ limits are indicated by the shaded region around the horizontal grey lines showing the amplitude of the injected signal. The four dominant peaks in the lower trace indicate successful recovery of all four signals, at amplitudes that are consistent with the injected values, and whose scatter is consistent with the expected uncertainty.
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.
For a set of n planet signals, the net orbital velocity vector v orb can be modelled as the product of a set of coefficient pairs θ orb = {A 1 , B 1 , · · · , A n , B n } with an array of time-domain function pairs F = {cos ω 1 t j , sin ω 1 t j , · · · , cos ω n t j , sin ω n t j }, ω k being the orbital frequency of the kth planet: For simplicity we assume circular orbits here, though eccentric orbits could in principle be fitted with periodic signals including additional Fourier components. 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 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 θ orb that minimizes χ 2 = (P ⊥ · δv T ) · −1 · (P ⊥ · δv), where δv ≡ v obs − v obs − F · θ orb .
Defining v ⊥ = P ⊥ · v obs and F ⊥ = P ⊥ · F, the goodness of fit is quantified by which is minimized by solving for θ orb : The log likelihood of the data given the model is Here we use the simplifying assumption that the RV measurements are uncorrelated. The covariance matrix is diagonal and its log determinant is l j =1 ln jj . The diagonal elements jj = Var(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 −1 . The sample mean and standard deviation (σ/ √ 4) of the signal amplitudes are 0.433 ± 0.005 m s −1 , 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 −1 semiamplitude 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 v obs at the same four periods, we obtain semiamplitudes 0.466, 0.560, 0.481, and 0.129 ms −1 , whose sample mean and standard deviation are 0.409 ± 0.083 m s −1 . Fig. 13 shows clearly the improvement in fidelity of the recovered amplitudes when the optimal shape model is applied. Table 1. Frequencies, periods, and semi-amplitudes of the strongest signals in the periodograms of raw and shape-corrected apparent velocities. The first two columns give the periods and semiamplitudes of the four injected signals. Columns 3-5 give the periods, semi-amplitudes, and uncertainties of all peaks with amplitudes greater than 33 cm s −1 in the periodogram of residual velocities remaining after subtraction of the SCALPELS projection, as in a blind RV search. The final three columns give the same information, from simultaneous modelling of CCF shape changes and planetary motion, made with prior knowledge of the four injected periods, as described in (Section 4.2).

Injected
Velocities from residual CCF Velocities from simultaneous fit  Figure 12. As for Fig. 11, but for the case where signal separation is performed simultaneously with orbit fitting given prior knowledge of the orbital periods. The middle (green) periodogram shows the difference between the observed and shape-driven velocities. The fourth (red) trace shows the periodogram of the fitted model of the four orbital signals. The bottom (magenta) trace shows the residuals after subtraction of both the shape-driven and orbital RV models. The scatter in the amplitudes of the four dominant peaks in the middle trace is consistent with the expected uncertainty, and their mean amplitude is unbiased relative to the injected values.
The formally propagated 0.066 m s −1 uncertainties in the semiamplitudes are ∼ 1.6σ v / √ N obs for N obs = 853 if we adopt a singlemeasurement precision of σ v 1.25 m s −1 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 Figure 13. The recovered amplitudes of the injected signals are shown in green for the simultaneous SCALPELS fit with k max = 13, and in blue for simultaneous sinusoidal fits to the uncorrected original RV data (k max = 0) at the injected periods. The grey line and band show the original signal level and the formal 1σ uncertainty on the recovered amplitude. The result demonstrates clearly the improvement in consistency and fidelity in the recovered signal amplitudes when SCALPELS is used. 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 timedomain 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.

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.

Planets injected into solar observations
Our algorithm recovered the RV signals of injected planets with semi-amplitudes of just 40 cm s −1 with an SNR of 6 based on a time series of 853 CCF measurements.
The semi-amplitude uncertainty of 6.6 cm s −1 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 ∼ 15 per cent precision for RV semiamplitudes as small as K = 40 cm s −1 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 semiamplitudes 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 −1 , 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.

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.

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 −1 signal with a period of 225 d, given an instrumental whitenoise uncertainty of 80 cm s −1 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 −1 and P = 300 d. Our results, using real data but an idealized observing pattern, give similar cause for optimism.

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 V ⊥ from 1.25 to 1.08 m s −1 . 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.

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 −1 amplitudes may be superposed on much larger signals of giant planetary companions at longer orbital periods. If the giant's orbit is wellcharacterized, 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.

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).

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.

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 −1 (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 v 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. U A (t) and 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.

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.

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).

AC K N OW L E D G E M E N T S
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.

DATA AVA I L A B I L I T Y
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.

A P P E N D I X A : C A L C U L AT I O N O F P RO F I L E D E R I VAT I V E S
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).
The numerical derivatives are calculated with a simple differencing scheme: Here h is the velocity increment per CCF sampling interval in velocity. The current HARPS-N DRS outputs CCFs on to a standard grid of l = 161 velocities with h = 0.25 km s −1 . Here C(v i , t j ), C (v i , t j ), and C (v i , t j ) are estimates of the CCF and its derivatives, which we describe below. 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.

A P P E N D I X B : TAY L O R -S E R I E S V E L O C I T Y M E A S U R E M E N T F RO M C C F
The RVs derived by the HARPS-N DRS are obtained by fitting a Gaussian profile to the activity-distorted CCF. Pure Doppler shifts simply displace the distorted profile in velocity. Following the formulation of Bouchy, Pepe & Queloz (2001) for line shifts that are small compared to the line width, we measure Doppler shifts using a first-order Taylor-series approximation to the instantaneous CCF profile shape at the time t j of the jth observation. For shifts much smaller than the 0.25 km s −1 velocity increment per CCF array element, the first derivative C (v i , t j ) of the instantaneous profile C suffices to measure the displacement of the mean-subtracted CCF The derivatives C (v i , t j ) are calculated from a reduced-rank representation of C(v i , t j ) as described in Appendix A. The covariance matrix is calculated as described in Appendix B1 below. The corresponding variances of the estimated RVs are given by .
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 = Cov(R) = R T · R/mestimates the systematic covariances between the columns of 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 −1 (Dumusque et al. 2020). The DRS delivers a default CCF that is also sampled in velocity increments of 0.82 km s −1 , 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.
The full covariance matrix can be computed via singular-value decomposition of R, as defined in equation (4): This approach gives results that are identical to a direct evaluation of the unweighted sample covariance matrix, with a maximum fractional deviation of 1 part in 1000. More importantly, it allows us to calculate a reduced-rank version of the covariance matrix, using only the leading k max principal components. With an appropriate choice of k max , a reduced-rank reconstruction of the mean-subtracted CCF residuals gives a representation of the large-scale temporal covariance pattern that can be written as 1 m R T kmax · R kmax = 1 m (P C,kmax · diag(S 2 C,kmax ) · P T C,kmax ).
This produces the low-pass filtered covariance matrix shown in the upper panel of Fig. B1. 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.
The profile of the ridge is seen clearly when the values of the covariance matrix are averaged along diagonals parallel to the leading diagonal, and plotted against velocity lag relative to the leading diagonal, as shown in the bottom panel of Fig. B1. Using a triangular fit to this average profile we obtain an average base half-width δv = 0.82 kms −1 . The base width of the ridge perpendicular to the diagonal is therefore about 2 CCF velocity increments, as expected from matching of the CCF sampling interval to the spectrograph resolution element and linear interpolation in the calculation of the CCF. We optimize k max to give a clear separation between the largescale column covariances and the covariances between neighbouring elements in each row arising from oversampling of photon noise. To Figure B1. Upper panel: reduced-rank (k max = 10) covariance matrix of the sequence of CCFs residuals, showing the large-scale covariance due to profile variability. Middle panel: Residual covariance matrix obtained by subtracting the reduced-rank representation from the full covariance matrix. Lower panel: Mean residual covariance along diagonals as a function of horizontal velocity offset from the leading diagonal. The diagonal ridge is triangular in crosssection, with base width 1.64 km s −1 , which is equivalent to two physical pixels on the spectrograph CCD. 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 −1 (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.
The row variances σ 2 j of R − R kmax reflect the SNR of the individual observations, so we use them as the scale factors A j = σ 2 j Figure B2. Apparent RVs derived from the CCF time series in the heliocentric frame via the Taylor approximation (equation B1) versus apparent RVs from the HARPS-N DRS after applying the barycentric to heliocentric correction. The mean values of both sets of velocities have been subtracted. The standard deviation of the Taylor-series velocities is 2.287 ms −1 . The two estimators of radial velocity are highly correlated and have slope near unity. The points are colour-coded with BJD-2400000.0. The RMS scatter in the differences between the two sets of measured velocities is 0.044 m s −1 .
of the triangular profile. This leads to the following model for the covariance matrix of the high-frequency noise of an observation made at time t j : This form of the covariance matrix is suitable for calculating the velocities and their precision. The resulting velocities differ from those obtained assuming spatially uncorrelated noise only at the 15 cm s −1 level. Their precision also depends, however, on the scatter introduced by profile-shape changes. When fitting an orbit to the velocities, it is therefore more appropriate to use velocity variances calculated with equation (B2) with a covariance matrix that also includes the reduced-rank model of the covariances arising from the time-varying shape of the CCF:

B2 Comparison with DRS velocities
We use this covariance matrix with equations (B1) and (B2) to estimate the RVs δ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 −1 .
The formal errors derived from equation (B2) are typically 0.11 m s −1 , though the data set includes a small number of days of sparse data yielding uncertainties greater than 0.2 m s −1 . 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 −1 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 = Cov(R kmax ) from equations (B3) and (B4). The resulting systematic error per observation is 0.82 ± 0.017 m s −1 , which is very close to the RMS amplitude of the scatter in Fig. 9.