# FIESTA II. Disentangling stellar and instrumental variability from
exoplanetary Doppler shifts in Fourier domain

J. Zhao Department of Astronomy & Astrophysics, The Pennsylvania State
University, 525 Davey Lab, University Park, PA 16802, USA  Exoplanetary
Science at UNSW, School of Physics, UNSW Sydney, NSW 2052, Australia  Eric B.
Ford Department of Astronomy & Astrophysics, The Pennsylvania State
University, 525 Davey Lab, University Park, PA 16802, USA  Center for
Exoplanets & Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State
University, University Park, PA 16802, USA  Center for Astrostatistics, 525
Davey Laboratory, The Pennsylvania State University, University Park, PA
16802, USA  Institute for Computational & Data Sciences, 525 Davey Laboratory,
The Pennsylvania State University, University Park, PA 16802, USA  C. G.
Tinney Exoplanetary Science at UNSW, School of Physics, UNSW Sydney, NSW 2052,
Australia  Australian Centre for Astrobiology, UNSW Sydney, NSW 2052,
Australia

###### Abstract

The radial velocity detection of exoplanets is challenged by stellar
spectroscopic variability that can mimic the presence of planets and by
instrumental instability that can further obscure the detection. Both stellar
and instrumental changes can distort the spectral line profiles and be
misinterpreted as apparent RV shifts.

We present an improved FourIEr phase SpecTrum Analysis (FIESTA a.k.a.
$\mathit{\phi}$ESTA) to disentangle apparent velocity shifts due to a line
deformation from a true Doppler shift. $\mathit{\phi}$ESTA projects stellar
spectrum’s cross correlation function onto a truncated set of Fourier basis
functions. Using the amplitude and phase information from each Fourier mode,
we can trace the line variability at different CCF width scales to robustly
identify and mitigate multiple sources of RV contamination. For example, in
our study of the 3 years of HARPS-N solar data, $\mathit{\phi}$ESTA reveals
the solar rotational effect, the long-term trend due to solar magnetic cycle,
instrumental instability and apparent solar rotation rate changes. Applying a
multiple linear regression model on $\mathit{\phi}$ESTA metrics, we reduce the
weighted root-mean-square noise from 1.89 m/s to 0.98 m/s. In addition, we
observe a $\sim$3 days lag in the $\mathit{\phi}$ESTA metrics, similar to the
findings from previous studies on BIS and FWHM.

Exoplanet detection methods (489), Radial velocity (1332), Fast Fourier
transform (1958), Stellar activity (1580), Astronomy data analysis (1858)

††software:  astropy (Astropy Collaboration et al., 2013, 2018), SOAP 2.0
(Dumusque et al., 2014), wpca (https://github.com/jakevdp/wpca), george
(Ambikasaran et al., 2015).

##  1 Introduction

In radial velocity exoplanet surveys, the Doppler shifts that result from the
gravitational interaction of planets and their host star are measured via
high-resolution spectroscopy to detect and characterise the masses and orbits
of planets. However, the spectral signature of a Doppler shift can be mimicked
by deformation of the host star’s spectral line profiles.

As the field pushes towards improved Doppler precision, developing robust and
powerful approaches for characterizing line profile variations, whether they
be due to intrinsic stellar variability or instrumental effects, is becoming
increasingly critical. Early planet-hunting spectrographs went to great
lengths to provide a precise and stable wavelength calibration, but did not
stabilise the instrument, leading to significant instrumental line spread
function variations. Recently, a new generation of spectrographs (e.g.,
ESPRESSO (Pepe et al., 2014), NEID (Schwab et al., 2019), the EXPRES
Spectrometer (Petersburg et al., 2020)) have been designed to be much more
intrinsically stable, in an effort to enable extremely precise radial velocity
(EPRV) measurements. These instrument specifications aim to reach 10-30 cm/s
radial velocity precision, with the long-term goal of characterizing
potentially Earth-like planets around Sun-like stars (Fischer et al., 2016).
To reach this goal, these instruments incorporate multiple environmental
control strategies, in an effort to stabilise the line-spread function and
minimise instrument-induced variability. Equally, or perhaps even more
importantly, this improved instrumental stability makes it feasible to use
line shape variations as a diagnostic to recognise and mitigate the effects of
stellar variability.

Astronomers have been exploring a variety of strategies to recognise and
mitigate the effects of stellar variability, ranging from looking for
correlations with traditional activity indicators such as the equivalent width
of H$\alpha$ (Herbst & Miller, 1989), the line profile full width at half-
maximum (FWHM, Queloz et al., 2009), the line profile’s bisector (Queloz et
al., 2001), and data-driven and machine learning methods (Pearson et al.,
2018; de Beurs et al., 2020, e.g.).

Zhao & Tinney (2020) proposed $\mathit{\phi}$ESTA for characterizing a signal
deformations and signal shifts. In the context of radial velocity detection of
exoplanets, the signal is usually the cross-correlation function (CCF) of the
stellar spectrum with a template spectrum or synthetic mask. By combining
information from many spectral lines, the CCF has significantly higher signal-
to-noise ratio (SNR) than individual spectral lines. In this study, we
validate and test $\mathit{\phi}$ESTA using CCFs based on a weighted binary
mask (Baranne et al., 1996; Pepe et al., 2002).

$\mathit{\phi}$ESTA decomposes the CCF into the orthogonal Fourier basis
functions and calculates the shift for each basis function. A pure Doppler
shift results in all the Fourier basis functions being shifted by the same
$RV_{0}$. In contrast, a signal varying its shape results in different basis
functions being shifted by different amounts. Zhao & Tinney (2020) used the
averaged effect of shifts in the lower and higher frequency ranges (denoted as
$RV_{\text{FT,L}}$ and $RV_{\text{FT,H}}$ respectively) as a summary statistic
to characterise the CCF behaviour.

In this paper, we provide a more rigorous derivation of the
$\mathit{\phi}$ESTA methodology, and make updates to ease the interpretation
of $\mathit{\phi}$ESTA outputs. For example, we abandon zero-padding (adding
zeros at the end of the signal to increase sampling) as was used in the
previous implementation of Fourier transform, since this acted as
interpolating the power spectrum and the phase spectrum but did not add useful
information. We also adopt the discrete Fourier transform (DFT), so that the
output sampling size in the velocity frequency domain matches the input in the
velocity domain, and thus for $N$ velocity grids in the CCF as input, there
will be $N$ outputs in the corresponding “velocity frequencies”.

The paper is organised in the following manner. We provide an intuitive way to
understand the maths behind the $\mathit{\phi}$ESTA method and discuss its
implementation in Section 2. We discuss the practical consideration when
applying $\mathit{\phi}$ESTA to the parametrisation of CCFs and analysing
timeseries in Section 3. In Section 4, we demonstrate the applications of
$\mathit{\phi}$ESTA on the SOAP 2.0 simulated solar spectra, with varying
latitudes of solar spots and plages, as well as the simulated solar
observation timeseries. Then it is followed by the applications of
$\mathit{\phi}$ESTA on the 3-years HARPS-N high-resolution solar observations
in Section 5. A detailed comparison of $\mathit{\phi}$ESTA models and similar
models using more traditional activity indicators (FWHM and BIS) is carried
out in Section 6. Section 7 is the summary and discussions where we also
compare $\mathit{\phi}$ESTA with previous works and discuss future research
opportunities.

Appendix A defines the coefficient of determination $R^{2}$ and the adjusted
$R^{2}$. Then in Appendix B we discuss the noise propagation in the Fourier
transform and explore under what conditions the distributions of
$\mathit{\phi}$ESTA amplitudes and phases are nearly Gaussian. Lastly,
Appendix C demonstrates that zero-padding (implemented in the earlier version
of $\mathit{\phi}$ESTA ) only changes the sampling of the Fourier transform
output and so we do not adopt it anymore.

##  2 Method

###  2.1 $\mathit{\phi}$ESTA in a nutshell

$\mathit{\phi}$ESTA is a Fourier decomposition of the CCF, whose flux
$CCF(v_{n})$ $(n=0,1,\dots,N-1)$ are evaluated at a discrete grid of
velocities, where the velocity grid $\\{v_{n}\\}$ is equally spaced and shared
across observations.

$CCF(v_{n})=\frac{1}{N}\sum_{k=0}^{N-1}\widehat{CCF}(\xi_{k})e^{i\frac{2\pi}{N}nk}\qquad
n=0,1,\dots,N-1.$ (1)

It is known as the inverse discrete Fourier transform (inverse DFT) that
decomposes $CCF(v_{n})$ into a linear combination of the orthogonal basis
functions $e^{i\frac{2\pi}{N}nk}$, weighted by the complex coefficient
$\widehat{CCF}(\xi_{k})$. Each coefficient $\widehat{CCF}(\xi_{k})$ is derived
from the DFT of $CCF(v_{n})$

$\widehat{CCF}(\xi_{k})=\sum_{n=0}^{N-1}CCF(v_{n})e^{-i\frac{2\pi}{N}nk}\qquad
k=0,1,\dots,N-1.$ (2)

The line shape information stored in $CCF(v_{n})$ is lossless in the
$\widehat{CCF}(\xi_{k})$ if all $N$ terms are retained, as both can be
converted via DFT and inverse DFT.

We define, for each mode $k$ in the Fourier domain representation, the
amplitude

$A_{k}=|\widehat{CCF}(\xi_{k})|,$ (3)

the phase

$\phi_{k}=\arctan{\frac{\widehat{CCF}(\xi_{k})_{\text{Im}}}{\widehat{CCF}(\xi_{k})_{\text{Re}}}}$
(4)

and the velocity frequency

$\xi_{k}=k/L$ (5)

where $L$ the length of the CCF velocity grid and $1/L$ is the unit velocity
frequency. Note that in this paper, the original domain is labeled in velocity
$(v)$ and the Fourier transformed domain is labeled in “velocity frequency”
($\xi$). Because the inputs $CCF(v_{n})$ are real functions and the DFT output
is Hermitian-symmetric, only the first $N_{0}=N/2+1$ or $(N+1)/2$ (depending
on whether $N$ is even or odd) modes are needed for CCF parametrisation and
the rest modes contain redundant information.

Replacing $\frac{nk}{N}$ by $\xi_{k}v_{n}$, we can rewrite the CCF in Eq. 1 in
velocity $v_{n}$ as opposed to index $n$

$CCF(v_{n})=\frac{1}{N}\sum_{k=0}^{N-1}\widehat{CCF}(\xi_{k})e^{i2\pi\xi_{k}v_{n}}\qquad
n=0,1,\dots,N-1.$ (6)

$\mathit{\phi}$ESTA interprets the CCF shapes, as parametrised by the
amplitudes $A_{k}$ and the phases $\phi_{k}$ at corresponding velocity
frequencies $\xi_{k}$ with $k=0,1,\dots,N_{0}-1$ in the Fourier domain. In
practice, we will keep only the leading terms which provide a dimensionally
reduced representation of the CCF. The decision of how many terms to retain
will depend on the spectral resolution, pixel sampling of the line spread
function and SNR. We discuss the choice for terms needed in Section 3.

We use the python function numpy.fft.rfft to compute the DFT of the real input
$CCF(v_{n})$ with the Fast Fourier Transform algorithm (Cooley & Tukey, 1965;
Press et al., 2007).

###  2.2 $\mathit{\phi}$ESTA applied to a pure line shift

A pure radial velocity shift $RV$ due to orbiting exoplanets results in a bulk
shift of the CCF as shown below.

$\displaystyle CCF(v_{n}-RV)$
$\displaystyle=\frac{1}{N}\sum_{k=0}^{N-1}\widehat{CCF}(\xi_{k})\cdot
e^{i2\pi\xi_{k}(v_{n}-RV)}$
$\displaystyle=\frac{1}{N}\sum_{k=0}^{N-1}\widehat{CCF}(\xi_{k})\cdot
e^{i(2\pi\xi_{k}v_{n}-2\pi\xi_{k}RV)}$ (7)

Comparing Eq. 7 with Eq. 6, the $RV$ shift results in a phase shift
$\Delta\phi_{k}=-2\pi\xi_{k}RV$, which is consistent with the differential
phase spectrum discussed in the continuous domain (Zhao & Tinney, 2020). While
the phase shift depends on the velocity frequency $\xi_{k}$, the ratio
$\Delta\phi_{k}/\xi_{k}$ is invariant across velocity frequencies and
proportional to the radial velocity shift.

$RV=-\frac{\Delta\phi_{k}}{2\pi\xi_{k}}.$ (8)

The amplitudes $A_{k}$ are also invariant for a pure line shift.

###  2.3 $\mathit{\phi}$ESTA applied to a perturbed line shape

Substituting $\widehat{CCF}(\xi_{k})$ by $A_{k}\cdot e^{i\phi_{k}}$ in Eq. 6,
we have the following form

$CCF(v_{n})=\frac{1}{N}\sum_{k=0}^{N-1}A_{k}\cdot e^{i\phi_{k}}\cdot
e^{i2\pi\xi_{k}v_{n}}.$ (9)

The shape of a deformed CCF can be characterised by changes in both the
amplitudes $A_{k}$ and the phases $\phi_{k}$, and thus

$\displaystyle CCF_{\text{deformed}}(v_{n})$
$\displaystyle=\frac{1}{N}\sum_{k=0}^{N-1}(A_{k}+\Delta A_{k})\cdot
e^{i(\phi_{k}+\Delta\phi_{k})}\cdot e^{i2\pi\xi_{k}v_{n}}$
$\displaystyle=\frac{1}{N}\sum_{k=0}^{N-1}(A_{k}+\Delta A_{k})\cdot
e^{i\phi_{k}}\cdot
e^{i2\pi\xi_{k}\big{(}v_{n}+\frac{\Delta\phi_{k}}{2\pi\xi_{k}}\big{)}}$ (10)

Comparing Eq. 10 with Eq. 9, the line deformations can be interpreted as an
apparent (but spurious) radial velocity shift for each velocity frequency $k$.

$RV_{\text{FT},k}=-\frac{\Delta\phi_{k}}{2\pi\xi_{k}}.$ (11)

In addition to the phase-derived $RV_{\text{FT},k}$, amplitude changes can
also be used to trace a perturbed line shape.

###  2.4 Measuring $\Delta\phi_{k}$ and $\Delta RV_{k}$

Measuring the phase shift $\Delta\phi_{k}$ requires specifying a reference
CCF. In this manuscript, the weighted mean CCF serves as our phase reference.
Throughout the paper (unless otherwise specified), we use the inverse variance
weighting to calculate the weighted averages (e.g., mean CCF, daily binned RV,
weighted RMS), where the variance is taken as the RV uncertainty squared. Once
we calculate $\Delta\phi_{k}$ between an observed and the reference CCF, Eq.
11 provides a group of relative radial velocity shifts $RV_{\text{FT},k}$
between the two CCFs for each velocity frequency $\xi_{k}$. Similar to
$RV_{\text{FT,L}}$ and $RV_{\text{FT,H}}$ proposed in Zhao & Tinney (2020)
that represent the RV shifts in the lower and higher frequencies,
$RV_{\text{FT},k}\,(k=0,1,\dots,N_{0}-1)$ will show the same response to a
bulk line shift but different responses to line shape changes. Therefore,

$\Delta RV_{k}\equiv RV_{\text{FT},k}-RV_{\text{apparent}}\qquad
k=0,1,\dots,N_{0}-1$ (12)

is used to parametrise the line profile deformation, where
$RV_{\text{apparent}}$ is a single apparent RV shift derived from the full
CCF. $RV_{\text{apparent}}$ may be computed using an independent RV
measurement algorithm or as a weighted mean of $RV_{\text{FT},k}$’s that
account for varying measurement precision (e.g. Section 4.2). In this paper,
we use $RV_{\text{Gaussian}}$, which is the centre of a Gaussian fit to the
CCF extending to $5\sigma$ in either direction from centre.

###  2.5 Summary

The effects of a pure shift (Eq. 8) and the effects of line shape changes (Eq.
11) are fundamentally different, even though the expressions are similar. The
$RV$ in Eq. 8 is the same for all velocity frequencies $\xi_{k}$, whereas the
$RV_{\text{FT},k}$ in Eq. 11 changes with the Fourier mode $k$. We treat Eq. 8
as a special case of Eq. 11 where $RV_{\text{FT},k}$ is a constant for all
$k\,(k=0,1,\dots,N_{0}-1)$. As a result, we can now use a group of
$RV_{\text{FT},k}\,(k=0,1,\dots,N_{0}-1)$ to parametrise the CCF behaviour (be
it a bulk shift, a deformation, or both). While a single measurement of the
apparent radial velocity shift per observation can not distinguish between the
effects of orbiting exoplanets and stellar variability, the set of
$RV_{\text{FT},k}$’s can be used to recognise the presence of line shape
changes. We explore that possibility using simulated and actual observations
in Section 4 and Section 5.

Changes in the Fourier amplitudes ($A_{k}$) also indicate line variability.
However, CCF shape changes in the vertical direction that may not induce
apparent RV variation can also impact $A_{k}$. In this paper, as we intend to
model the line variability induced RV, we focus on the analysis of
$RV_{\text{FT},k}$ and $\Delta RV_{k}$, which are directly related to the
spurious RV (see Section 4.2 for further discussions).

##  3 Useful Fourier modes in $\mathit{\phi}$ESTA

Noise in the CCF propagates into the uncertainties in the $\mathit{\phi}$ESTA
amplitudes and phase (and thus radial velocity shifts calculated from Eq. 11).
We need to understand the effects of noise in order to interpret the results
of the $\mathit{\phi}$ESTA analysis. In principle, using all $N$ Fourier modes
(or $N_{0}$ Fourier modes for real inputs as in our case) can reconstruct the
CCF without losing information. In practice, measurement noise and/or the
spectrograph resolution may limit the utility of higher frequencies. In this
section, we consider each of the factors that may limit the number of useful
terms in the $\mathit{\phi}$ESTA reconstruction. We summarise five aspects to
consider when it comes to making use of the amplitudes $A_{k}$ and phases
$\phi_{k}$ information in the presence of associated uncertainties in the
Fourier domain as below.

  1. 1.

How many Fourier modes are necessary to accurately and adequately reconstruct
the typical CCF shape?

  2. 2.

What is the distribution of errors in $A_{k}$ and $\phi_{k}$ due to photon
noise?

  3. 3.

From which $k$ does the SNR of $A_{k}$ and $\phi_{k}$ become so low that the
estimates of individual $A_{k}$ and $\phi_{k}$ are no longer useful?

  4. 4.

From which $k$ is the variability in the $A_{k}$ and $\phi_{k}$ timeseries
dominated by the measurement uncertainty of $A_{k}$’s and $\phi_{k}$’s?

  5. 5.

Which terms of the $\mathit{\phi}$ESTA decomposition correspond to velocity
scales greater than the width of the spectrograph line spread function (LSF)?

They are discussed in the following subsections respectively.

###  3.1 Reconstructing the CCF

First, we consider how many Fourier modes are needed for $\mathit{\phi}$ESTA
to accurately parametrise a typical CCF without overfitting. The number of
modes needed depends on the shape of the CCF and the modelling error, which is
defined as the residual root mean square (RMS) between the observed CCF and
the $\mathit{\phi}$ESTA reconstruction of the CCF using only the first $k$
terms. Too few terms would result in an inaccurate representation of the CCF
due to a large residual, but too many terms would result in a modelling error
below the photon noise, indicating overfitting.

For example, Fig. 1 shows the first 5 $\mathit{\phi}$ESTA basis functions (or
6 including the zeroth mode, which is the mean function of the CCF) used to
fit an arbitrary but typical HARPS-N solar spectrum CCF. For details of the
solar spectrum, refer to Dumusque et al. (2021) and Section 5. Using up to
$k=5$ $\mathit{\phi}$ESTA basis functions reaches the residual or modelling
error of $6.7\times 10^{-4}$, corresponding to a SNR111It is approximately
calculated as the median flux of the CCF between -15 and 15 km/s (signal)
divided by the modelling error (noise). = $1.4\times 10^{3}$ for the CCF it is
able to represent. The HARPS-N solar CCFs have a high median SNR $\sim$ 10,000
for individual CCFs and SNR $\sim$ 46,000 for the daily binned CCFs. In order
for the modelling error to be less than the photon noise level, one would need
to use all the available Fourier modes, which is $k_{\textrm{max,photon}}=20$.
Note that $k_{\textrm{max,photon}}$ is recommended for extracting all usable
information of the CCF up to the photon noise limit, but any $k\leq
k_{\textrm{max,photon}}$ provides valuable information.

Figure 1:  (a) The $\mathit{\phi}$ESTA basis functions used to fit a HARPS-N
solar spectrum CCF up to $k=5$ Fourier modes. (b) The residual of the CCF
after fitting the first $k$ Fourier modes. (c) The modelling precision (i.e.
residual RMS) as a function of the number of modes up to $k$ used to
parametrise the CCF. The vertical axis on the right stands for the SNR to
which the CCF can be reconstructed within the photon noise uncertainties.
Using 20 Fourier modes will fully recover the CCF (i.e., zero modelling noise;
not shown). The modelling precision as a function of $k$ is for illustration
only and it may vary if a different CCF profile is used, e.g., from a
different star or a different instrument LSF.

###  3.2 Noise distribution of amplitudes and phases

According to Eq. 2, $\widehat{CCF}(\xi_{k})_{\mathrm{Re}}$ and
$\widehat{CCF}(\xi_{k})_{\mathrm{Im}}$ can be obtained by a linear projection
of $CCF(v_{n})$ onto a set of Fourier basis vectors. Therefore, Gaussian noise
in $CCF(v_{n})$ translates into Gaussian errors for
$\widehat{CCF}(\xi_{k})_{\mathrm{Re}}$ and
$\widehat{CCF}(\xi_{k})_{\mathrm{Im}}$. In fact, as derived in Appendix B,
their $1\sigma$ uncertainties are

$\sigma_{\widehat{CCF}(\xi_{k})_{\text{RE}}}=\sigma_{\widehat{CCF}(\xi_{k})_{\text{IM}}}=\sqrt{\frac{1}{2}\sum_{n=0}^{N-1}\sigma^{2}_{CCF(v_{n})}}\coloneqq\widehat{\sigma}$
(13)

where $\sigma_{CCF(v_{n})}$ is the uncertainty of $CCF(v_{n})$.

In contrast, the amplitude, phase and the resulting $RV_{\text{FT},k}$ are
obtained by non-linear transformations (i.e., square and square root
operations on two random variable in Eq. 3 and the ratio and arctan operations
in Eq. 4). Therefore, the uncertainties in $A_{k}$’s, $\phi_{k}$’s and
$RV_{\text{FT},k}$’s are not strictly Gaussian. The amplitudes tend to be
skewed with a lower boundary at zero and the phases tend to present a narrower
core and broader tails (Fig. 15). These effects can be substantial when the
input CCFs data are noisy. We define the noise level (NL) as

$\text{NL}=\frac{\widehat{\sigma}}{A_{k}}.$ (14)

Based on numerical simulations, we have verified that noise in measurement of
amplitudes and phases follow nearly the normal distribution when NL $\leq$ 0.2
(Appendix B). Under this low NL regime, we can approximately derive the
analytical expression for the uncertainties of amplitudes and phases

$\sigma_{A_{k}}=\widehat{\sigma}\quad\text{and}\quad\sigma_{\phi_{k}}=\frac{\widehat{\sigma}}{A_{k}}.$
(15)

In general, the higher the mode $k$ becomes, the smaller the amplitudes and
thus the larger the NL. We define a cut-off frequency index based on normality
$k_{\text{max,normal}}$ and $\xi_{k_{\text{max,normal}}}$ as the last
frequency that satisfies our criterion. By limiting our analysis to $k\leq
k_{\max,\text{normal}}$, we ensure that noise in $A_{k}$ and $\phi_{k}$
closely follows a normal distribution.

###  3.3 Individual $A_{k}$ and $\phi_{k}$ SNR

For an individual observation, we compute the SNR for $A_{k}$’s and
$\phi_{k}$’s as $A_{k}/\sigma_{A_{k}}$ and $1/\sigma_{\phi_{k}}$ 222As the
phase $\phi_{k}$ ranges between $-\pi$ and $\pi$, a larger $|\phi_{k}|$ does
not necessarily mean the signal is stronger. We therefore use the unity
instead of $|\phi_{k}|$ in the numerator to calculate the phase SNR. using the
uncertainties from that observation. The SNR for $A_{k}$’s and $\phi_{k}$’s
are thus identical according to Eq. 15. As $A_{k}$ decreases rapidly while
$\sigma_{A_{k}}$ is nearly a constant (see Eq. 13 and 15) with increasing $k$,
the SNR for $A_{k}$’s and $\phi_{k}$’s decreases for higher $k$. Therefore, we
limit our analysis to $k\leq k_{\max,\textrm{individual}}$, where
$k_{\max,\textrm{individual}}$ is the largest $k$ for which the median SNR of
$A_{k}$’s or $\phi_{k}$’s is at least 2.

When choosing which $k$ to analyse in a timeseries, we use the median
$\sigma_{A_{k}}$ and median $\sigma_{\phi_{k}}$ instead of the uncertainty for
each observation, so as to ensure consistency across observations and to avoid
discarding information in $k$-th timeseries that would result in only a small
fraction of low SNR observations not passing our minimum SNR threshold.

###  3.4 Timeseries SNR

When applying $\mathit{\phi}$ESTA to a timeseries of CCFs, we aim to
characterise the deviations of ${A_{k}}$ and $\phi_{k}$ from their time-
averages. We define a $\mathrm{SNR}_{\mathrm{timeseries},k}$ as the ratio of
the RMS deviation of each $A_{k}$ and $\phi_{k}$ from their mean value over a
timeseries to the median $\sigma_{A_{k}}$ and $\sigma_{\phi_{k}}$.

It is useful to limit our analysis to the $A_{k}$’s and $\phi_{k}$’s that have
measurement uncertainties less than the extent of their variations with time.
In this paper, we require a minimum SNR of 2 to ensure high-quality timeseries
are obtained, i.e., the standard deviation of an $A_{k}$ or $\phi_{k}$
timeseries is at least twice as large as the median $\sigma_{A_{k}}$ or median
$\sigma_{\phi_{k}}$ in the timeseries. We define $k_{\max,\text{timeseries}}$
as the maximum $k$ for which this criterion is satisfied.

###  3.5 Instrumental resolution

Finally, we consider $k_{\max,\mathrm{inst}}$, the maximum $k$ for which we
expect changes in the CCF could be dominated by the target star, rather than
the instrument or detector. Due to the instrumental broadening of the spectral
lines, higher Fourier modes that captures changes below the LSF width will not
be useful in parametrising stellar variability. For example, the instrumental
resolution of HARPS in velocity space is 2.5 km/s. For a CCF with a velocity
ranging from $-15$ to 15 km/s, Fourier modes higher than
$k=L/(1/\xi_{k})=30/2.5=12$ (Eq. 5) capture variations on a scale finer than
the instrumental resolution. In this example, we limit the $A_{k}$’s and
$\phi_{k}$’s analysed to have $k\leq k_{\max,\mathrm{inst}}=12$. As shown in
Section 5, the solar rationally induced variability are barely seen in modes
with $k>12$.

###  3.6 Summary notes

For illustrative purposes, we use the three years of HARPS-N solar CCFs
(details in Section 5) to give readers an idea of what $k$ is recommended
based on the criteria discussed above.

$k=\min(k_{\textrm{max,photon}},k_{\text{max,normal}},k_{\max,\textrm{individual}},k_{\max,\text{timeseries}},k_{\max,\mathrm{inst}})=\min(20,20,20,11,12)=11$

will ensure that the $A_{k}$’s and $\phi_{k}$’s being analysed can be
associated with an unbiased $1\sigma$ uncertainty and have sufficient SNR to
characterise the CCF and the timeseries variations as well as focusing on
revealing information that could be due to stellar variations (as opposed to
instrumental or detector variations). While $\mathit{\phi}$ESTA is capable of
including higher $k$’s, its default choice of thresholds (and those used in
this study) are designed to be conservative in ensuring that the measurement
uncertainties are very nearly Gaussian.

##  4 Validation of $\mathit{\phi}$ESTA on simulated solar observations

We use the SOAP 2.0 (Spot Oscillation And Planet 2.0, Dumusque et al., 2014)
simulator to generate CCFs of solar spectra affected by spots or plages that
are based on the HARPS observations of the Sun. Unless specified in this
paper, we set the spot at $30^{\circ}$ latitude with the size of
$0.1R_{\sun}$; the stellar configurations are based on the Sun and chosen as
the default values of SOAP 2.0, e.g. $T_{\text{star}}=5778~{}K$,
$T_{\text{spot-photosphere}}=663~{}K$.

###  4.1 Line shift vs line deformation

First we generate the CCFs of solar spectra affected by a single spot in one
rotation. In the presence of stellar variability, each $RV_{\text{FT},k}$
differs for each CCF, because each mode $k$ traces the line deformation at a
different length-scale, with $k=1$ tracing variability at the full range of
the CCF, $k=2$ tracing variability at 1/2 the range, etc. Therefore,
$RV_{\text{FT},k}$ provides a quantitative measurement of the radial velocity
shift and the line deformation.

To illustrate the changes in the Fourier amplitudes and phases of the CCFs for
a line deformation as opposed to a line shift, we take one snapshot at solar
rotation phase 0.51, when the spot is slightly off solar disk centre and the
resulting CCF is deformed with an apparent radial velocity measured at 1.1
m/s. For comparison, we shift the undeformed CCF by the same 1.1 m/s and feed
them into $\mathit{\phi}$ESTA. When the CCF is decomposed into the Fourier
basis functions, the amplitudes and RV shifts remain the same for all modes
for a pure line shift, while amplitudes and RV shifts vary for a line
deformation (Fig. 2).

Figure 2: Amplitudes and phases at a glance for two CCFs, one caused by a pure
shift (blue dashed line) and the other caused by a spot at $30^{\circ}$
latitude with the size of $0.1R_{\sun}$ at rotation phase 0.51 (red solid
line). $\Delta\textrm{CCF}$, $\Delta A_{k}$, $\Delta\phi_{k}$ and $\Delta
RV_{k}$ are relative measurements to the template CCF which is neither shifted
or deformed. Even though both RVs are comparable (1.1 m/s),
$\mathit{\phi}$ESTA tells them apart (one is intrinsic and the other is
spurious) and provide a quantitative measurement of the phase and amplitude
changes. Note that it is impossible to replicate a line profile with a shift
below the sampling spacing, and so even the pure shift of a CCF simulated by
interpolation will end up in a slightly different shape and result in unequal
RV shifts in higher frequencies when the amplitudes are smaller.

###  4.2 Weighted mean of $RV_{\text{FT},k}$

We calculate the weighted mean of $RV_{\text{FT},k}$ for the SOAP simulated
CCFs over the course of one rotation, with the weights being the corresponding
amplitudes $A_{k}$. We find that the weighted means are consistent with the
radial velocity shifts measured by fitting a Gaussian function to the CCF
(Fig. 3). Therefore, the apparent radial velocity $RV_{\text{Gaussian}}$ can
be treated as the overall effect of the radial velocity shifts of individual
Fourier basis functions (i.e. $RV_{\text{FT},k}$). $RV_{\text{FT},k}$, on the
other hand, can be treated as a multi-dimensional measurement of the RV shift
at different CCF length scales.

Figure 3: Left panel: timeseries of the first 5 $RV_{\text{FT},k}$ over the
course of solar rotation for SOAP 2.0 simulations of a single solar spot of
$0.1R_{\sun}$ at $30^{\circ}$ latitude. Only rotation phases between 0.2 and
0.8 are shown. Right panel: the weighted mean of $RV_{\text{FT},k}$
overplotted with $RV_{\text{Gaussian}}$. It shows while individual
$RV_{\text{FT},k}$ traces the line deformation differently, the overall effect
of $RV_{\text{FT},k}$ is consistent with $RV_{\text{Gaussian}}$, indicating
the set of $RV_{\text{FT},k}$ is a multi-dimension measurement of the RV shift
that has an overall effect equivalent to $RV_{\text{Gaussian}}$.

###  4.3 Spots vs plages and their $\mathit{\phi}$ESTA correlations

We further explore how $\mathit{\phi}$ESTA behaves under the influence of a
solar spot and plage in various latitudes (from 10 to 80 degrees) in one
stellar rotation. Earlier studies using the previous version of
$\mathit{\phi}$ESTA can be referred to Zhao (2019).

As we can tell from Fig. 4, higher latitudes tend to have minor effects on
$A_{k}$ ($\Delta A_{k}$ is shown using a solar CCF without spot or plage as a
reference) and $\Delta RV_{k}$. This is because both the projected area of
spots and plages in high latitudes and their rotational velocity along the
line of sight are smaller. On the other hand, changes in $A_{k}$ and $\Delta
RV_{k}$ per latitude towards the equator are not as prominent as in high
latitudes, and this is because a latitude change near the equator does not
change the the projected area of spots and plages along the line of sight and
the rotation speed as much as in higher latitudes.

Figure 4: $\Delta A_{k}$ and $\Delta RV_{k}$ in one solar rotation for spots
(left two panels) and plages (right two panels), simulated in latitudes from
10 to 80 degrees.

Since changes in both $A_{k}$ and $\Delta RV_{k}$ are parametrisations for the
amount of spectral line variability, a larger deformation, which is likely to
lead to larger apparent radial velocity shifts, is expected to show larger
changes in $A_{k}$ and larger $\Delta RV_{k}$ as well. However, this is not
always the case as the correlations presented in Fig. 5. Overall, larger
apparent RV tend to be only associated with larger $\Delta RV_{k}$ for spots
in lower Fourier modes (e.g. $k=1$ and 2). For plages and higher Fourier
modes, their linear correlations are not observed. On the flip side, the
different patterns of $\Delta RV_{k}$ \- $RV_{\text{Gaussian}}$ plots and
their stronger dependence on latitude will be valuable for figuring out the
latitude information of the spots and plages.

Figure 5: $\Delta RV_{k}$ and the apparent RV correlation for a simulated
solar spot (top 5 panels) and a solar plage (bottom 5 panels) in one solar
rotation. The latitudes are classified by the rainbow colours. Linear
correlations between $\Delta RV_{k}$ and $RV_{\text{Gaussian}}$ are only
prominent for spots in lower Fourier modes ($k=1$ and 2). Spot and plage
configurations are the same as the ones in Fig. 4.

###  4.4 Simulated continuous solar observations

In the next step, we analysed SOAP 2.0 solar simulations based on the known
spot properties considering a similar activity level to the Sun, including
spot temperature, size distribution and lifespan (Gilbertson et al., 2020a).
The simulated data set was sampled twice per day and has 730 spectra covering
a whole year. It was originally used for timeseries analysis (Gilbertson et
al., 2020b), we nevertheless use the simulated solar CCFs as individual
observations and seek correlations between $\Delta RV_{k}$ and the apparent RV
due to solar spots.

We present the RV and $\Delta RV_{k}$ timeseries of the simulated data set in
Fig. 6. The results are consistent with the discussions in Section 4.3 that
$\Delta RV_{k}$ shows strong linear correlations with the apparent RV
($RV_{\text{SOAP}}$) for $k=1$ and $2$, as the (Gilbertson et al., 2020b)
solar simulations are deformed only by spots but not plages. It indicates the
$\mathit{\phi}$ESTA lower Fourier modes may be used for linear RV correction
for spot-dominated stars.

Figure 6: RV and $\Delta RV_{k}$ timeseries of SOAP 2.0 simulations of the Sun
with evolving spots over a year. Units are in m/s. The lower Fourier modes
($k=1$ and 2) show strong linear correlations with $RV_{\text{SOAP}}$, as
quantified by the adjusted coefficient of determination $\bar{R}^{2}$
(Appendix A). The corresponding periodograms are on the right, with red dots
labelling the prominent periodicities and vertical lines labelling
periodicities near the solar rotation and its harmonics at 26.6, 14.6 and 8.0
days.

##  5 $\mathit{\phi}$ESTA on HARPS-N 3-years solar observations

###  5.1 Data

As an example to show how $\mathit{\phi}$ESTA can analyse real world line
profile variability, we apply $\mathit{\phi}$ESTA to the high-cadence HARPS-N
spectroscopic observations of the Sun-as-a-star during 2015 to 2018 (Dumusque
et al., 2015; Phillips et al., 2016). The data are publicly accessible on
https://dace.unige.ch/sun/, including 34550 observations. These public data
sets were already pre-selected to minimise non-astrophysical effects on
observed CCF (Collier Cameron et al., 2019). They were taken in good observing
condition (i.e. not significantly affected by clouds), the differential
extinction for each observation less than 0.1 m/s, and the data reduction
software (DRS) reports good quality results. The same selection criteria were
also used in de Beurs et al. (2020) and Dumusque et al. (2021). The HARPS-N
CCFs and RVs used in Dumusque et al. (2021), Collier Cameron et al. (2021) and
in this paper (daily binned as in first row of Fig. 7) were processed by the
new HARPS-N DRS as opposed to those used in de Beurs et al. (2020) processed
by an older version of the HARPS-N DRS. As a result, the $RV_{\textrm{HARPS}}$
presented in this paper reveals a long-term RV drift induced by the solar
magnetic cycle and contributes to a larger scatter of apparent RVs (both
before and after applying mitigation procedures) than the ones in de Beurs et
al. (2020).

We work on the inverted, normalised daily weighted averaged CCF and the the
daily weighted RVs, therefore reducing the number of original observations
over the three years from 34550 to 567. Each exposure is designed to integrate
over the solar p-mode oscillation timescale of $\sim$ 5 min (Strassmeier, K.
G. et al., 2018; Chaplin et al., 2019). The daily weighted average CCFs
further average out both the p-mode solar oscillation effects and much of
solar granulation which is significant on timescales from minutes for granules
(Vázquez Ramió, H. et al., 2005) to hours for mesogranules (Matloch, L. et
al., 2009). A $\mathit{\phi}$ESTA study on these intra-day oscillation and
granulation effects will be carried out in the near future.

The periodogram of the daily weighted HARPS-N RVs (first row of Fig. 7 or 8)
shows several prominent periodicities in the daily RV timeseries, with
$\sim$28.5 days being identified as the solar rotation period. The other
periods will be discussed in the following sections. The unbinned RV
periodogram (not shown) is almost identical to the daily binned RV periodogram
above 1 day.

###  5.2 $\mathit{\phi}$ESTA timeseries

In the left-hand panels of Fig. 7 and 8, we show the timeseries for each
$\Delta RV_{k}$ and $A_{k}$ for 20 Fourier modes. Different modes characterise
the CCF deformation on different scales, with $k=1$ being the longest mode
(refer to Fig. 1). As $k$ increases, the scale of study decreases as $1/k$ and
the amplitude decays approximately as a Gaussian function of $k$ (Fig. 2).
Therefore, the signal-to-noise ratio for measurements of $\Delta RV_{k}$ and
$A_{k}$ become dominated by photon noise for large $k$. The middle panels are
the correlations between each row’s indicator and $RV_{\text{HARPS}}$. The
right-hand panels of Fig. 7 and 8 show the periodogram for each $\Delta
RV_{k}$ and $A_{k}$.

Based on the $\Delta RV_{k}$ timeseries and the timescales observed in their
periodograms, we’ve identified five mechanisms contributing to the variations
in the CCF profile:

  1. 1.

Solar rotationally modulated CCF variations: The $\sim 28$ day solar rotation
period and its harmonics can be seen in both $\Delta RV_{k}$ and $A_{k}$
periodograms (e.g. $k$ from 3 to 8). The relative significance of the peaks
often differs between the periodograms of $\Delta RV_{k}$ and $A_{k}$, even
for the same $k$. Rotationally modulated variations are more obvious in
$\Delta RV_{k}$ periodograms.

  2. 2.

Solar magnetic cycle: The long-term RV drift over the three years of
observations is correlated with $\Delta RV_{k}$ and $A_{k}$ up to $k=4$
($\bar{R}^{2}\sim 0.5$).

  3. 3.

CCD detector warm-ups: Instrumental features can be seen most clearly as the
sudden jumps of $\Delta RV_{1}$ and $\Delta RV_{2}$ that coincide with times
of interventions to the CCD detector (Dumusque et al., 2021) . While not
strictly periodic, these this effect likely leads to the periodogram peaks
near $\sim 210$ d and $\sim 285$ d.

  4. 4.

Changes in the apparent solar angular velocity: The Earth-Sun distance changes
both due to Earth’s eccentricity and as a result of perturbations from Jupiter
(also discussed in Collier Cameron et al., 2019). The periodicity at around
400 days is predominant in the higher Fourier modes (e.g. $k\geqslant 9$) and
nearly matches the time between conjunctions with Jupiter.

  5. 5.

Changes in the apparent solar angular velocity as a result of the Earth’s
orbit being slightly eccentric changing the distance between Earth and the Sun
(Collier Cameron et al., 2019). This periodicity at half a year is apparent in
many of the higher Fourier modes, but can be overwhelmed by the CCD detector
warm-ups and/or harmonics of the timescale between conjunctions of Jupiter.

Figure 7: Left: RV inferred from HARPS-N solar observations (top row) and
$\Delta RV_{k}$ for the first 20 Fourier modes. $RV_{\textrm{HARPS}}$ and
$\Delta RV_{k}$ timeseries have units of m/s. Middle: The correlations between
each $\Delta RV_{k}$ and $RV_{\textrm{HARPS}}$ described by the adjusted
coefficient of determination $\bar{R}^{2}$. Right: The periodograms of the
corresponding timeseries with prominent periods labelled in red dots. The
vertical bands correspond to periods at 13.4, 28.5, 365/2, 200-220, 285 and
400 days as discussed in Section 5.2. Figure 8: Same as Fig. 7, but for
$A_{k}$ instead of $\Delta RV_{k}$.

While both $\Delta RV_{k}$ and the change in $A_{k}$ trace the CCF
variability, $A_{k}$ seems to capture more periodicities. This is expected
since line depth changes or artifacts in the estimate for the continuum around
the line may affect $A_{k}$’s, but not $\Delta RV_{k}$’s.

###  5.3 Principal component analysis

Some of the nearby Fourier modes present similar or correlated behaviour as
seen in the timeseries and their periodograms, especially for $\Delta RV_{1}$
and $\Delta RV_{2}$, which appear anti-correlated with each other. It is
possible that certain features of a line deformation can take place across
several Fourier modes and thus cause correlated effects. Therefore, we can
model a large fraction of the CCF variance with a smaller number of
coefficients by invoking principal component analysis (PCA) to perform
dimensional reduction on the $\Delta RV_{k}$’s (or $A_{k}$’s). This has dual
benefits of increasing the SNR in the retained PC scores and potentially ease
interpretation of a lower dimensional timeseries.

In an ideal world with ultra-high SNR, we could make use of all the 20 Fourier
modes for the CCF parametrisation. In reality, higher Fourier modes are more
strongly affected by measurement uncertainties; in addition, they do not
contribute as much information due to rapidly decaying amplitudes (Section 3).
As mentioned in Section 3.6, it is advised to use no more than 11 Fourier
modes for the HARPS-N solar data analysis. If we were to focus on the first 3
principal components built from the first $k_{\textrm{max}}$ modes, we find
that increasing $k_{\textrm{max}}$ from 7 to 11 hardly changes the first 2 PC
score timeseries, and increasing $k_{\textrm{max}}$ from 9 to 11 hardly
changes the $3^{\textrm{rd}}$ PC score timeseries. Therefore, we choose
$k_{\textrm{max}}=9$ to parametrise the CCF variability. When performing the
analysis in Sections 5.5 and 5.6, we also find that building the first 3
principal components from the first 9 Fourier modes results in the smallest
residual (compared with other $k_{\textrm{max}}$), though the difference of
the residual WRMS (i.e. weighted RMS) is less than 5% among choices of
$k_{\textrm{max}}$ ranging from 7 to 11.

While most PCA packages are designed to deal with equally weighted data by
default, the HARPS-N observations have associated errorbars. The weighted PCA
is aimed to tackle this problem by solving the weighted covariance matrix of
the data (Delchambre, 2014). We used the wpca python package for its
implementation. The weighted PCA on $\Delta RV_{k}$ (Fig. 9) reduces the
feature dimension from 9 to 3 while $\sim$75% of variance in the data are
preserved. The residual due to projecting data into a lower dimension accounts
for the rest $\sim$25% of the variance. We can tell the aforementioned
periodicities in Section 5.2 are preserved in the first 3 PC scores. Note that
we normalise the $\Delta RV_{k}$ for each $k$ before applying weighted PCA to
prevent the higher modes from dominating the principal components.

Figure 9: Left: PCA for $\Delta RV_{k}$ up to $k_{\textrm{max}}=9$. The first
3 principal components explain 75% of the variance in $\Delta RV_{k}$
$(k=1,2,\dots,9)$. Middle and the right columns are the same correlations
plots and periodograms plots as in Fig. 7 except for these principal
components.

###  5.4 Separating out long-term and short-term variations

PCA extracts independent features in high dimensional data while keeping most
of the variation information within the first few principal components,
however, the process is not physically driven and thus we should not expected
that each mechanism that contributes to CCF variations will be neatly
contained in an individual principal component. Therefore, we apply a low-pass
filter to the PC scores to separate rotationally induced effects (periods
$<100$ days) from long-term effects (periods $>100$ days). Specifically, we
compute L-PC${}_{i}(t)$, the mean of Gaussian process conditioned on each of
the PC score timeseries (PCi’s) separately using a Matérn 5/2 kernel (Genton,
2001; Minasny & McBratney, 2005, e.g.) with a length-scale of $r=100$ days,
which captures the variations above $r$ and smooths out those below $r$ (Fig.
10). We then subtract L-PCi from the corresponding principal component scores
to obtain S-PCi, which probes the short-term variability below $\sim$100 days
and are found to be dominated by the solar rotationally modulation.

Figure 10: Separating the short-term (S-PC) and long-term (L-PC) variability
of the 3 principal components in Fig. 9. The corresponding periodograms are on
the right.

###  5.5 Multiple linear regression modelling

Once PCA has reduced the number of features required to accurately parametrise
the CCF variation, we want to investigate whether these PC scores are useful
for predicting the behaviour of the spurious RV measurements. The simplest
approach would be a multiple linear regression model

$y\equiv RV_{\text{HARPS}}=\mathbf{X}\beta+\epsilon$ (16)

where $\mathbf{X}$ is the regressor matrix that traces the CCF variability,
$\beta$ is the coefficient vector, $\epsilon$ is the error vector and $y\equiv
RV_{\text{HARPS}}$ is the response vector.

We examine two models for which $\mathbf{X}$ consists of various features
derived from the CCFs parametrisation, Model 1 using the first three principal
component scores (PCi’s computed from $\Delta RV_{k}$, $k=1,2...,9$, Fig. 9)
and Model 2 taking a step further by using their short-term components
$\textrm{S-PC}_{i}(i=1,2,3)$ and long-term components
$\textrm{L-PC}_{i}(i=1,2,3)$ as regressors (Fig. 10). As models naturally fit
the data better with more features or parameters, to avoid over-fitting, we
perform Lasso regression (i.e., added a L1 regularisation penalty to the
weighted sum of the squared residuals, Tibshirani, 1996; Hastie, 2015) and
minimise the following

$\text{loss
function}=(y-\mathbf{X}\beta)^{T}\mathbf{W}(y-\mathbf{X}\beta)+\lambda\|\beta\|_{1}$
(17)

where $\mathbf{W}$ is the weight matrix, $\lambda$ is the tuning parameter and
$\|\beta\|_{1}$ is the 1-norm of the vector $\beta$ defined as
$\sum_{i=1}^{n}|\beta_{i}|$. $\lambda=0$ implies no penalty and Eq. 17 is
reduced to the normal linear regression optimisation. The larger $\lambda$
becomes, the larger the penalty term and the more coefficients are driven to
zero. To find out the optimal $\lambda$, we adopt a cross-validation approach
where the data are randomly scrambled and split into 5 folds. We train our
model as we minimise the loss function in Eq. 17 using 4 folds of the RV data
and feature matrix, and then obtain a residual WRMS for the validation data.
This process is repeated 100 times, each with a randomised split and $\lambda$
is tested from $10^{-3}$ to 1. Then we locate $\lambda=\hat{\lambda}_{1}$
corresponding to the smallest WRMS in the testing set. Usually a more
regularised and simpler model with $\lambda>\hat{\lambda}_{1}$ is preferred.

With the multiple linear regression model regularised by the loss function
above (Eq. 17), we reduce the uncorrected WRMS of $RV_{\textrm{HARPS}}$ from
2.00 m/s to 1.36 m/s for using the first 3 PC scores (Fig. 9) and further down
to 1.14 m/s with additionally separating out the long-term and short-term
variations (Fig. 10), corresponding to 32% and 43% reduction in the WRMS, or
53% and 67% reduction in the weighted variance, as calculated by
$\bar{R}^{2}$. Among the total reduction in the weighted variance, we
determine the importance of each feature by calculating their contributed
variance $(\beta_{i}\sigma_{i})^{2}$. The information is summarised in Table
1’s rows corresponding to Models 1 and 2.

We estimate that solar rotation-induced RVs contributes just $\sim 7\%$ of the
variance (based on sum of the variance percentages due to
$\textrm{S-PC}_{i}(i=1,2,3)$) or $\sim 0.4$ m/s in terms of RMS without
considering measurement errors or instrumental instability. In contrast,
$\textrm{L-PC}_{1}$ (which traces the solar magnetic cycle and the CCD
detector warm-ups) contributes over 80% of the variance reduction. Most of the
remaining $\sim$10% of the RV variation is believed to be caused by changes in
the apparent solar rotation rate due to Jupiter and Earth’s eccentricity
($\textrm{L-PC}_{2}$).

Model | Label | $\beta_{i}$ | $\sigma_{i}$ | $(\beta_{i}\sigma_{i})^{2}$ | Data WRMS | Model WRMS | Residual WRMS | $\bar{R}^{2}$
---|---|---|---|---|---|---|---|---
1 | $\textrm{PC}_{1}$ | 0.73 | 1.96 | 2.06 (99.9%) | 2.00 | 1.44 | $1.36~{}(-32\%)$ | 0.532
$\textrm{PC}_{2}$ | -0.00 | 1.25 | 0.00 (0.0%)
$\textrm{PC}_{3}$ | 0.04 | 1.06 | 0.00 (0.1%)
2 | $\textrm{S-PC}_{1}$ | 0.51 | 0.71 | 0.13 (4.8%) | 2.00 | 1.55 | $1.14~{}(-43\%)$ | 0.669
$\textrm{S-PC}_{2}$ | 0.23 | 1.05 | 0.06 (2.1%)
$\textrm{S-PC}_{3}$ | -0.03 | 0.81 | 0.00 (0.0%)
$\textrm{L-PC}_{1}$ | 0.83 | 1.81 | 2.24 (80.7%)
$\textrm{L-PC}_{2}$ | -0.87 | 0.63 | 0.30 (10.8%)
$\textrm{L-PC}_{3}$ | 0.35 | 0.62 | 0.05 (1.7%)
3 | (see Fig. 11) | 1.89 | 1.53 | $0.94~{}(-50\%)$ | 0.733
4 | FWHM | 0.29 | 0.95 | 0.08 (4.0%) | 2.00 | 1.43 | $1.34~{}(-33\%)$ | 0.552
BIS | 1.35 | 1.00 | 1.79 (96.0%)
5 | S-FWHM | 0.70 | 0.70 | 0.24 (10.4%) | 2.00 | 1.54 | $1.18~{}(-41\%)$ | 0.651
S-BIS | 0.00 | 0.48 | 0.00 (0.0%)
L-FWHM | 0.16 | 0.59 | 0.01 (0.4%)
L-BIS | 1.66 | 0.86 | 2.06 (89.2%)
6 | (see Fig. 13) | 1.89 | 1.46 | $1.04~{}(-45\%)$ | 0.680
Table 1: Comparisons of 6 models fitting the HARPS-N solar spurious RV. The
difference among the models are the regressors used in the multiple linear
regression model. Note that the adjusted coefficient of determination
$\bar{R}^{2}$ accounts for the the number of free parameters in a model and so
$\bar{R}^{2}$ can indicate the performance of a model regardless of different
number of parameters (Appendix A).
Model 1: using PC scores derived from $\Delta RV_{k}$ ($k_{\textrm{max}}=9$).
Model 2: in additional to Model 1, using the solar rotation components
$\textrm{S-PC}_{i}$ and the long-term variation components $\textrm{L-PC}_{i}$
derived from the PC scores are separated out.
Model 3: in additional to Model 2, introducing time lags for the solar
rotation components $\textrm{S-PC}_{i}$.
Models 4 - 6: repeating for Models 1-3 but for FWHM and BIS.

###  5.6 Multiple linear regression modelling with time lags

Collier Cameron et al. (2019) explored the correlation between the spurious RV
signal in HARPS-N observations and two traditional CCF shape indicators – the
Bisector Inverse Slope (BIS) and the full-width half-maximum (FWHM). They
found that the apparent RV variations led the BIS by 3 days and the FWHM by 1
day. We want to investigate if introducing a lag between $RV_{\textrm{HARPS}}$
and the principal components of $\Delta RV_{k}$ can improve our predictions
for the spurious RV signal. As the lag between two timeseries can be poorly
constrained when the variation timescale are too long compared to the time
lag, we decide to fix a zero lag for the long-term variations
$\textrm{L-PC}_{i}$ and only introduce a lag parameter in integer days for the
short-term variations $\textrm{S-PC}_{i}$.

In order to implement multiple linear regression with time lags, we must
address two challenges. First, observations are not available for every
consecutive day. Most gaps are less than a few days whilst a large one starts
from 2017-11-11 (BJD$-2400000=58068.5$) due to the HARPS-N Fabry-Pérot upgrade
(Dumusque et al., 2021). Second, the daily binned timestamps are not exactly
evenly spaced, creating a slight misalignment between timeseries when they are
shifted by integer number of days. To overcome these challenges, we linearly
interpolate the indicators timeseries to the $RV_{\textrm{HARPS}}$ timestamps
and their offsets in days. The beginning and the end of the two trunks of
$RV_{\textrm{HARPS}}$ separated by the Fabry-Pérot upgrade gap are not used so
that we can avoid the artifacts of extrapolating the indicators timeseries
outside the $RV_{\textrm{HARPS}}$ timestamps when dealing with lags. We also
note that over 90% of the daily binned observations are within 2 hours of the
day, so using the $RV_{\textrm{HARPS}}$ directly as the response vector is
sufficient to study the time lags in days.

The lagged multiple linear regression introduces additional fitting parameters
– one coefficient for each lag for each of the $\textrm{S-PC}_{i}$ and
$\textrm{L-PC}_{i}$ ($i=1,2,3$). We tested different maximum lags allowed and
found that the fitted coefficients remain mostly stable once the model
includes at least a maximum 5-day lag. Therefore, we present the results using
a maximum 5-day lag. As for the regularisation parameter $\lambda$, we note
the residual WRMS reaches minimum at $\hat{\lambda}_{1}=0.028$. According to
“the one standard error rule” (Breiman et al., 1984; Hastie et al., 2009), a
more regularised model with $\hat{\lambda}_{2}=0.080$ at one standard error of
$\textrm{WRMS}(\hat{\lambda}_{1})$ is usually recommended (Fig. 11). However,
considering we would like to compare all the models using $\bar{R}^{2}$, it
would be easier if the same $\lambda$ is chosen. Therefore, we choose the
middle ground $\lambda=0.05$, which results in a more parsimonious Model 3
than if $\lambda=\hat{\lambda}_{1}$ were chosen and meanwhile $\lambda=0.05$
does not over regularise Models 1 and 2333For Models 1 and 2, their residual
WRMS curves remain mostly flat until they start rising at $\lambda\sim 0.01$,
indicating hardly any overfitting and so regularisation is not needed.
However, in order to fairly compare different models, we choose $\lambda=0.05$
for all the models.. With this, we obtain a residual WRMS = 0.94 m/s, a 50%
reduction in the WRMS or 73% reduction in the weighted variance.

Figure 11: Left: we determine the tuning parameter $\lambda$ in the Lasso
regression model (Eq. 17) as follows. Let $\lambda=\hat{\lambda}_{1}$
correspond to the smallest WRMS in the testing set (red solid line; with the
red dashed indicating the $1\sigma$ uncertainty). Generally, a more
regularised and simpler model with a larger $\lambda$ up to
$\hat{\lambda}_{2}$ is preferred, where
$\text{WRMS}(\hat{\lambda}_{2})=\text{WRMS}(\hat{\lambda}_{1})+\sigma_{\text{WRMS}(\hat{\lambda}_{1})}$.
We choose $\lambda=0.05$ in order to be consistent with other models. Right:
the variance in percentage for each $\textrm{S-PC}_{i}$ and
$\textrm{L-PC}_{i}$ ($i=1,2,3$) in the regularised multiple linear regression
fitting of $RV_{\textrm{HARPS}}$. The blanks indicate zero coefficients (thus
zero variance) due to the regularisation with $\lambda=0.05$, while the 0.0’s
are the result of rounding from a non-zero, but small enough variance. A
positive lag means the proxy indicator lags behind the RV variation, i.e., the
RV variation leads the proxy indicator. The model returns $\bar{R}^{2}=0.733$
and a residual WRMS of 0.94 m/s.

In Fig. 11, we observe that $\sim 7\%$ of the variance (or 0.51 m/s in terms
of RMS) comes from the solar rotationally induced RV variability (i.e., sum of
the variance in $\textrm{S-PC}_{i}(i=1,2,3)$ for all lags), consistent with
the $\sim 7\%$ in Model 2 with no lags. Over 80% of variability is dominated
by terms that show trends of the solar magnetic cycle plus the CCD detector
warm-ups. In addition, the lag information of $\textrm{S-PC}_{i}(i=1,2,3)$
tracing solar rotation is spread out between $-2$ days and 5 days, with the
single most prominent lag for 3 days by summing up the variance of
$\textrm{S-PC}_{i}(i=1,2,3)$ by day, accounting for $\sim 2\%$ of the RV
variance. Individually, the weighted averaged lag for
$\textrm{S-PC}_{i}(i=1,2,3)$ are 2.2, 1.7 and 3.7 days, while Collier Cameron
et al. (2019) noted lags of 3 and 1 days between the BIS and FWHM and the
spurious RV. Our finding that a weighted average of the indicators at multiple
lags can provide a more effective predictor of the spurious RV signals is
complementary to results from Collier Cameron et al. (2019). We hypothesise
the different mechanisms that result in RV variation as the Sun rotates may
show different lags, and $\textrm{S-PC}_{i}(i=1,2,3)$, BIS and FWHM detect
part or a combination of these mechanisms.

##  6 The alternative: BIS and FWHM

To investigate if using $\mathit{\phi}$ESTA metrics improves the modelling of
line profile variability, we compare the results of §5.5 and 5.6 with a
similar analysis applied to two traditional CCF variability indicators, the
FWHM and BIS. We repeat the process of HARPS-N RV modelling except for
substituting the PC scores of $\Delta RV_{k}$ by the daily binned FWHM and BIS
(Fig. 12). The originally unbinned FWHM and BIS data were provided on
https://dace.unige.ch/sun/. The FWHM and BIS timeseries are standardised (mean
= 0 and variance = 1) so that the linear regression coefficients are not
penalised inappropriately in the Lasso regression because of their different
units or because of one having intrinsically smaller variance but requiring a
larger coefficient to compensate.

Figure 12: Same as Fig. 10 but for FWHM and BIS.

Comparing Fig. 10 and Fig. 12, $\textrm{L-PC}_{1}$ derived from
$\mathit{\phi}$ESTA appears to behave similarly to L-BIS, the long-term
variation of the BIS timeseries. As Fig. 13 shows, this long-term components
contribute to $\sim 91\%$ of variance in modelling $RV_{\textrm{HARPS}}$. The
other $\sim 9\%$ (equivalently 0.55 m/s RMS) comes from the short-term
variations which we attribute to be primarily solar rotationally modulated
variability. In terms of performance, using FWHM and BIS together results in a
slightly worse $\bar{R}^{2}=0.680$ and a slightly larger residual weighed RMS
of 1.04 m/s. This is likely due to the fact that the combination of FWHM and
BIS is not as effective at describing the CCF as the combination of
$\textrm{S-PC}_{i}$ and $\textrm{L-PC}_{i}$ derived from $\Delta RV_{k}$’s.

Figure 13: Same as Fig. 11 but for FWHM and BIS. The regularised model with
$\lambda=0.05$ returns $\bar{R}^{2}=0.680$ and a residual WRMS of 1.04 m/s.

We summarise 3 models using FWHM and BIS for predicting the HARPS-N solar
spurious RV in Table 1 (Models 4-6). In order to compare how the models are
similar to each other, we compare the residual correlations between Model 3
and Model 6 specifically. Their residuals are well correlated with each other,
indicating that the two models agree well on each other in predicting
$RV_{\textrm{HARPS}}$. The characteristic $1\sigma$ width of the residuals
(the half-width of the central 68.2% credible interval) is 0.77 m/s (Fig. 14),
indicating the overall difference of the two model predictions. We suggest
using this approach for evaluating the consistency between two or more models
in the planet detection or stellar variability modelling, such as model
comparisons in the EXPRES Stellar-Signals Project (Zhao et al., 2022).

Figure 14: Left: the correlation and histograms of residuals (i.e. the
difference between the observations and the model predictions) using the Model
3 and Model 6 predictions show that the two models are consistent. Right: the
difference between Model 3 and Model 6 predictions has a $1\sigma$ width of
0.77 m/s in the distribution.

##  7 Summary and discussions

###  7.1 Overview

We presented the improved Fourier phase spectrum analysis (FIESTA or
$\mathit{\phi}$ESTA) for disentangling intrinsic RV shifts from apparent RV
shifts caused by stellar variability and instrumental instability.
$\mathit{\phi}$ESTA provides summary statistics for each spectra in a
timeseries based on projecting the CCF onto a truncated version of the Fourier
basis functions. This approach is motivated by the shift-invariant properties
of the Fourier basis:

  1. 1.

A shift in the signal (a spectral line profile or a CCF in the context of
radial velocity exoplanet detection) does not change the amplitudes of each
Fourier basis function:

$A_{k}(t_{i})=A_{k}(t_{0})\qquad\textrm{for each $k$}$ (18)

  2. 2.

A true Doppler shift results in the same shift for each of the basis
functions:

$RV_{\text{FT},k=1}=RV_{\text{FT},k=2}=\dots=RV_{\textrm{shift}}$ (19)

and thus

$\Delta RV_{k=1}=\Delta RV_{k=2}=\dots=0.$ (20)

As a result, we can use changes in the amplitude $A_{k}$ and RVs measured as
$RV_{\text{FT},k}$ or $\Delta RV_{k}$ for each Fourier mode $k$ to parametrise
the variations of spectral line shapes. We focus on the analysis of $\Delta
RV_{k}$ because it is directly related to the spurious RVs due to stellar
and/or instrumental variability.

Specifically, $\mathit{\phi}$ESTA decomposes the spectral line profile CCF
into the Fourier modes (Eq. 6), i.e. the Fourier basis functions truncated to
have finite support (Fig. 1, left panel). The amplitude ($A_{k}$) and the
phase-derived RV shift ($RV_{\text{FT},k}$) for each Fourier mode $k$ fully
parametrise the CCF through discrete Fourier transform. The suite of
$RV_{\text{FT},k}$ is a $N$-dimensional measurement of the shared intrinsic
planetary RV shift of the CCF plus any stellar or instrumental variability
measured at mode $k$ (Section 4.2). The changes with time of the amplitudes
and RV shifts for each Fourier mode (denoted as $A_{k}$ and $\Delta RV_{k}$)
is used to trace the CCF variability.

As discussed in Section 3, higher Fourier modes generally come with larger
uncertainties. Practical consideration to pick the number of modes for the
analysis include (1) the precision of reconstructing the CCF, (2) the
uncertainties of $A_{k}$ and $RV_{\text{FT},k}$ for a single observation, (3)
the SNR of $A_{k}$ timeseries and $\Delta RV_{\text{FT},k}$ timeseries and (4)
the instrument resolution.

We demonstrated the use of $\mathit{\phi}$ESTA in analysing the RV timeseries
from Section 4 to Section 6, covering SOAP 2.0 solar simulations and the
3-years of HARPS-N solar observations.

##### SOAP 2.0 simulations

We demonstrated how $\mathit{\phi}$ESTA distinguishes an intrinsic line shift
from an apparent line shift due to a line deformation (Section 4.1). We also
explored the response of $\Delta RV_{k}$ to simulated solar spots and plages
at different latitudes, as well as their correlations with the apparent RV
(Section 4.3). The simulated solar observations show high correlations between
the spot-induced spurious RVs and the lower Fourier modes (Section 4.4).

##### HARPS-N solar data

In Section 5, we applied $\mathit{\phi}$ESTA to 3 years of HARPS-N solar
observations (Fig. 7 and 8). As neighbouring Fourier modes may show high
correlations with each other, we are motivated to use PCA for dimension
reduction. With PCA, we extracted the most prominent 3 orthogonal features
from the first 9 $\Delta RV_{k}$ timeseries (Fig. 9). Next, we separated the
short-term variability dominated by solar rotation from the long-term
variability modelled by a Matérn 5/2 GP kernel with a length scale of 100 days
(Fig. 10). Feeding these features into multiple linear regression models to
fit $RV_{\textrm{HARPS}}$, we reduced the WRMS RV from 2.0 m/s to 1.14 m/s
($\bar{R}^{2}=0.67$) in a model that uses only the spectra at each epoch.
Incorporating $\mathit{\phi}$ESTA outputs a few days prior/after each
observations allowed us to further reduce the WRMS RV from 1.89 m/s to 0.94
m/s ($\bar{R}^{2}=0.73$). To avoid overfitting, Lasso regression and cross
validation have been implemented. Our model shows the solar rotationally
induced variability contributes only 7% to the total RV variation modeled by
$\mathit{\phi}$ESTA (Fig. 11). The long-term variability, which we identified
as sources from the solar magnetic cycle and instrumental instability, compose
of over 80% of the RV variations in $RV_{\textrm{HARPS}}$. The remaining less
than 10% RV variability come from the CCF changing resulted from the apparent
solar rotation rate changes due to Jupiter and Earth’s eccentricity.
Therefore, it is important for future EPRV exoplanet surveys to consider
modelling and correcting for instrumental instability and long-term activity
variations, in addition to modelling the rotationally-linked stellar
variability.

Furthermore, we compared the models using PC scores of $\Delta RV_{k}$ as
features (Section 5) and the ones using FWHM and BIS (Section 6) as summarised
in Table 1. The latter performs slightly worse with a residual WRMS of 1.18
m/s ($\bar{R}^{2}=0.65$) without lags and 1.04 m/s ($\bar{R}^{2}=0.68$) with
lags. The fact that FWHM and BIS do not perform as well may be because they do
not provide as complete CCF parametrisation as the $\mathit{\phi}$ESTA-derived
$\Delta RV_{k}$. However, both methods predict the solar rotational spurious
RMS RV $\sim 0.5$ m/s. As a measurement of model consistency, the RMS of the
difference between the two residuals using these two models is $\sim 0.77$ m/s
(Fig. 14). We encourage future studies to perform a similar consistency check
for multiple methods proposed to reduce the effects of stellar variability.

###  7.2 Comparison to previous works

####  7.2.1 Original $\mathit{\phi}$ESTA

The original $\mathit{\phi}$ESTA proposed by Zhao & Tinney (2020) extracted
the CCF variability into the lower and higher frequency ranges
$RV_{\text{FT,L}}$ and $RV_{\text{FT,H}}$. The current version of
$\mathit{\phi}$ESTA presented here decomposes the CCF into all the available
Fourier modes up to the CCF sampling limit, providing a comprehensive
parametrisation of the CCF using the amplitudes $A_{k}$ and the phase derived
RV shift $RV_{\text{FT},k}$ and the CCF variability using $A_{k}$ and $\Delta
RV_{k}$. This facilitates the study of multiple sources contributing to CCF
variability, as demonstrated in our analysis of the HARPS-N solar observations
in Section 5.

####  7.2.2 SCALPELS

Collier Cameron et al. (2021) introduced SCALPELS to study the CCF variability
using the autocorrelation function (ACF) of the cross-correlation function of
each spectrum with a mask. SCALPELS was able to reduce the RMS of daily
averaged HARPS-N solar spectra from 1.76 m/s to 1.25 m/s, a 29% reduction in
the RV scatter. In order to make a more direct comparison to this result, We
conducted a preliminary analysis of the same HARPS-N solar observations as
Collier Cameron et al. (2021). When applied to the 5 year timeseries,
$\mathit{\phi}$ESTA reduced the weighted daily RV RMS to 1.09 m/s using the
same approach as in Section 5.6. Since the Collier Cameron et al. (2021)
dataset provides only daily binned CCFs and does not include FWHM and BIS
measurements and it became available after much of our analysis had been
completed, this manuscript only focuses on the three year timeseries available
on https://dace.unige.ch/sun/ at the time writing.

In addition to comparing the amount of stellar variability removed, it is
interesting to compare the basis vectors used by the two studies. The
eigenvectors for the ACF in Collier Cameron et al. (2021) Fig 3 look very
similar to the Fourier basis functions used in $\mathit{\phi}$ESTA (Fig. 1).
It is exciting to see how the two completely different methods - the data
driven approach from Collier Cameron et al. (2021) and the theoretical
derivation from Zhao & Tinney (2020) and this paper end up with similar basis
functions (i.e. eigenvectors).

####  7.2.3 Machine learning

de Beurs et al. (2020) analysed 3 years HARPS-N solar observations and predict
the spurious RV due to stellar variability by applying either linear
regression or a dense neural network to the CCF. Using a neural network, they
reduced the RV scatter by $47\%$, very similar to our $48\%$ reduction in the
WRMS using the multiple linear regression with $\mathit{\phi}$ESTA indicators
and allowing for a lag (Section 5.6). A precise comparison of these two
results is not practical because the de Beurs et al. (2020) HARPS-N RVs were
based on an earlier version of the HARPS-N data reduction system, while ours
are from the newer version (see discussions in Section 5.1). Regardless, there
is room for further research in reducing the residual scatter by treating the
$\mathit{\phi}$ESTA inputs (either directly or after dimensional reduction via
PCA) with neural network methods that are able to learn non-linear
relationships between activity indicator proxies and RVs.

###  7.3 Opportunities for Future Research

####  7.3.1 PCA

We explored PCA for dimension reduction and used the PC scores as features to
characterise the CCF deformation. As PCA is an unsupervised learning
technique, we could not separate out the different mechanisms that caused the
CCF deformation. In order to improve interpretability, we computed a smoothed
version of the PC scores ($\textrm{L-PC}_{i}$’s) that was insensitive to
changes on the solar rotation timescale. Indeed, this resulted in the long-
term solar magnetic cycle and the CCD detector warm-up events becoming more
prominent in $\textrm{L-PC}_{1}$. Inspecting the difference between the PC
scores and smoothed PC scores allowed us to separate out the short-term
variability (e.g. the solar rotation). In future studies, we could explore
employing regularisation in PC scores so as to drive some basis vectors to
focus on variations on a timescale near the stellar rotation period and other
basis vectors to focus on CCF shape changes occurring on longer timescales
(e.g. instrumental variations or long-term stellar cycles).

####  7.3.2 Gaussian processes

We used activity indicator proxies (PC scores derived from $\Delta RV_{k}$;
FWHM and BIS) as regressors in an effort to model the spurious RV introduced
by solar variability. Although they depend only on each instantaneous
spectrum, we showed that incorporating temporal information by performing
regression on these activity indicator proxies and the lagged components
simultaneously resulted in improved prediction power, as evidenced by
$\bar{R}^{2}$. This implies Gaussian process (GP) regression, particularly,
multivariate GP models with a common latent process and its derivatives
approximating the effect of regressing on the lagged indicator proxies would
be a promising approach. (Rajpaul et al., 2015; Jones et al., 2020; Gilbertson
et al., 2020b, GLOM,. Although GP may not be needed for the highly sampled
solar data analysed in this manuscript, planet hunting in other stellar
targets with sparse observing cadence would benefit from GP analysis of both
the stellar spurious RVs and the activity indicator proxies (e.g. Langellier
et al., 2021; Faria et al., 2022). In fact, we have made initial attempts to
combine the PC scores derived from $\Delta RV_{k}$ and the multi-variate GP
known as GLOM for modelling the stellar variability for HD 101501, HD 34411,
HD 217014, and HD 10700 as part of the EXPRES Stellar-Signals Project (Zhao et
al., 2020). This approach allows $\mathit{\phi}$ESTA to be applied to stars
other than the Sun, for which gaps in the observations are much more common
(Zhao et al., 2022).

Most published applications of GP regression to radial velocity timeseries
have assumed stational GP kernels. Stationary kernels are not well suited for
modelling sudden changes in the $RV_{\textrm{HARPS}}$ due to detector warm-ups
or power failures. We find that these contribute a significant fraction of the
CCF deformation changes in the HARPS-N solar dataset. Future studies may
benefit from adopting more complex GP kernels, such as the sum of a stationary
GP kernel for modelling stellar variability and a non-stationary kernel with
parameters informed by prior knowledge about any instrumental changes.

####  7.3.3 CCF construction

In this paper, we demonstrate the $\mathit{\phi}$ESTA methodology using
HARPS-N solar observations and simulated SOAP data. Publicly available CCFs
for HARPS-N solar data and SOAP data are based on a weighted binary mask.
While CCFs based on a template spectrum can improve the robustness of RVs for
some stars (e.g., low SNR, broad absorption features in cool stars), we expect
a narrow, weighted binary mask to preserve more line shape information. While
it could be interesting to study the effects of different CCFs masks and
preprocessing (e.g., CCF continuum normalisation, velocity range) on
$\mathit{\phi}$ESTA outputs, that is beyond the scope of this paper. Of
course, one should use the CCFs constructed using the same mask and
preprocessing to correctly account for other changes due to stellar
variability and instrumental instability.

####  7.3.4 Improving $\Delta RV_{k}$

$\Delta RV_{k}$, the difference between the RV shift at $k$-th Fourier mode
($RV_{\text{FT},k}$) and the apparent RV ($RV_{\mathrm{apparent}}$), traces
the amount of CCF variability at the frequency $\xi_{k}$ (Eq. 12). By
construction, a true Doppler shift ($RV_{\text{planets}}$) does not affect the
expected values of individual $\Delta RV_{k}$’s, but contributes to each of
the $RV_{\text{FT},k}$. For the HARPS-N solar spectra, there are no true
Doppler shift signals in the heliocentric RV timeseries and so
$RV_{\text{FT},k}$ alone would trace the CCF variability and could have been a
stronger indicator for variability. Nevertheless, for other stars, there will
inevitably be uncertainty about the potential presence of planetary
companions. That is why we choose to use $\Delta RV_{k}$ instead of
$RV_{\text{FT},k}$ even for the analysis for the HARPS-N solar observations as
a manner of consistency.

It may be possible to improve the statistical power and accuracy of $\Delta
RV_{k}$ as a series of CCF variability indicators by re-defining

$\widehat{\Delta RV_{k}}^{(l)}\equiv
RV_{\text{FT},k}-RV_{\text{planets}}^{(l)}\qquad k=0,1,\dots,N-1$ (21)

once we have a first guess of $RV_{\text{planets}}^{(1)}$, the RV signal of
putative planetary companion(s). Such first guess may be derived from jointly
fitting the planetary RV signal and intrinsic stellar variability with $\Delta
RV_{k}$. Then as an iterative approach, we keep improving the fit of
$RV_{\text{planets}}^{(l+1)}$ as we improve $\widehat{\Delta RV_{k}}^{(l)}$
from the previous round.

For its implementation, one could take a Bayesian approach and simultaneously
analyse the planetary signal and stellar variability. This would require
placing priors on both the planet properties and the $\mathit{\phi}$ESTA
indicators (e.g. $A_{k}$’s and $\Delta RV_{k}$’s). Multivariate Gaussian (GP)
regression (Section 7.3.2) provides a powerful framework for placing priors on
the $\mathit{\phi}$ESTA indicators. It is beyond the scope of this paper and
offers opportunities for future research.

##  8 acknowledgments

The authors gratefully acknowledge Suvrath Mahadevan, Jessi Cisewski Kehe,
Nikunj Sura and the Penn State team for their direct and indirect input to the
methodology, coding and writing of this manuscript. The analysis of this work
used the public three years of HARPS-N solar data hosted by University of
Geneva. We thank Xavier Dumusque and Zoe de Beurs for their opinions on the 3
years HARPS-N solar observations. We also thank Andrew Collier Cameron for the
productive discussions on SCALPELS and the 5 years HARPS-N solar observations.
J.Z. is grateful to Petra Brčić for the inspirational discussions on the
mathematical part of the $\mathit{\phi}$ESTA method over the course of
preparing the paper. J.Z. thanks Ludovic Delchambre and Jake Vanderplas for
clarifying the weighted PCA theories and implementations.

This research was supported by Heising-Simons Foundation Grant #2019-1177.
E.B.F. acknowledges the support of the Ambrose Monell Foundation and the
Institute for Advanced Study. This work was supported by a grant from the
Simons Foundation/SFARI (675601, E.B.F.). We acknowledge the Penn State Center
for Exoplanets and Habitable Worlds, which is supported by the Pennsylvania
State University and the Eberly College of Science.

We acknowledge the use of the following software for producing the results of
this research:

The Github repository for the analysis linked in the paper including the
$\mathit{\phi}$ESTA code is https://github.com/jinglinzhao/FIESTA-II, which is
also available in Zenodo (Zhao, 2022).

##  Appendix A The coefficient of determination $R^{2}$ and the adjusted
$R^{2}$

The coefficient of determination, also known as $R^{2}$, is defined as
$1-\frac{\textrm{wRSS}}{\textrm{wTSS}}$ where wRSS is the weighted sum of
squares of residuals and wTSS is the total weighted sum of squares (Heinisch,
1962). It is a measurement of the percentage of variance of the dependent
variable predicted by the independent variable(s), with $R^{2}=1$ indicating a
zero residual and the modelled values totally match the observations.

To avoid overestimating when extra explanatory variables are introduced into
the model, the adjusted $R^{2}$, denoted as $\bar{R}^{2}$, is used (Raju et
al., 1997). $\bar{R}^{2}$ is defined as $1-(1-R^{2})\frac{n-1}{n-p-1}$ where
$n$ is the sample size and $p$ is the number of explanatory variables.
$\bar{R}^{2}$ can thus be used to compare the performance of different models,
regardless of different explanatory variables being used. For the simple
linear regression, there is only a single explanatory variable and so $p=1$.

##  Appendix B Noise in amplitude and phases

For simplicity, we treat photon noise across the velocity grid as uncorrelated
and Gaussian (approximated for Poisson distribution because of high photon
counts), with the mean being the observed flux $CCF(v_{n})$ and the $1\sigma$
of the Gaussian distribution being the associated photon noise error counts at
each $v_{i}$ denote as $\sigma_{CCF(v_{n})}$. Then the uncertainty of
$\widehat{CCF}(\xi_{k})$ can also be calculated using DFT

$\sigma_{\widehat{CCF}(\xi_{k})}=\sum_{n=0}^{N-1}\sigma_{CCF(v_{n})}e^{-i\frac{2\pi}{N}nk}\qquad
k=0,1,\dots,N-1.$ (B1)

Since the DFT is a linear transform, a normally distributed $CCF(v_{n})$
results in a normally distributed $\widehat{CCF}(\xi_{k})$. The distribution
of $\widehat{CCF}(\xi_{k})$ can be approximated by a 2D multivariate normal
distribution in the Fourier domain. The marginalised
$\widehat{CCF}(\xi_{k})_{\text{RE}}$ and $\widehat{CCF}(\xi_{k})_{\text{IM}}$
also follow the normal distribution, with the variance

$\sigma^{2}_{\widehat{CCF}(\xi_{k})_{\text{RE}}}\approx\sigma^{2}_{\widehat{CCF}(\xi_{k})_{\text{IM}}}\approx\frac{1}{2}\sum_{n=0}^{N-1}\sigma^{2}_{CCF(v_{n})}.$
(B2)

Next we consider what distributions the amplitudes $A_{k}$ and $\phi_{k}$ as
calculated in Eq. 3 and 4 follow. We study under what condition $A_{k}$ and
$\phi_{k}$ are well-approximated by the normal distribution. Without loss of
generality, we choose $\widehat{CCF}(\xi_{k})=(1,0)$ and test in a
bootstrapping approach if the distribution of $A_{k}$ and $\phi_{k}$ can be
distinguished from a normal distribution when the noise-to-signal-ratio
$N/S=\sigma_{\widehat{CCF}(\xi_{k})_{\text{RE}}}=\sigma_{\widehat{CCF}(\xi_{k})_{\text{IM}}}=0.1,~{}0.2,~{},\dots,~{}1$
using the D’Agostino‘s K-squared (D’agostino et al., 1990) and the Shapiro-
Wilk normality tests (Shapiro & Wilk, 1965). We find that when $N/S\leq 0.2$,
the resulting $A_{k}$ and $\phi_{k}$ distributions are indistinguishable from
a normal distribution. Larger $N/S$ results in deviations from the normal
distribution (Fig. 15).

Figure 15: The distributions of the amplitudes and phases in the presence of
noise, where the signal (S) refers to $|\widehat{CCF}(\xi_{k})|=1$ and the
noise (N) refers to $\sigma_{\widehat{CCF}(\xi_{k})_{\text{RE}}}$ and
$\sigma_{\widehat{CCF}(\xi_{k})_{\text{IM}}}$. The black curve is a normal
distribution fit to the histogram (regardless whether the underlying
distribution is normal). N/S = 0.2 is the threshold when the distributions of
$A_{k}$ and $\phi_{k}$ start to deviate from normal distributions under the
normality test. For comparison, examples of N/S = 0.4 and 0.8 represent the
lower SNR regime, where the $A_{k}$ distributions tend to be skewed and the
$\phi_{k}$ distributions tend to be taller in the core and thinner in the wing
than a normal distribution.

In conclusion, we compute the uncertainties of
$\sigma_{\widehat{CCF}(\xi_{k})_{\text{RE}}}$ and
$\sigma_{\widehat{CCF}(\xi_{k})_{\text{IM}}}$, which can be approximated by
$\sqrt{\frac{1}{2}\sum_{n=0}^{N-1}\sigma^{2}_{CCF(v_{n})}}$, and compare it
with $|\widehat{CCF}(\xi_{k})|$. When the ratio is less than 0.2, we can be
confident that the distribution for the amplitudes and phases is well
described by a normal distribution.

Since we treat noise across CCF pixels uncorrelated, the resulting
uncertainties in the amplitudes and phases obtained by the bootstrapping are
likely to be overestimated compared to a noise model that accounts for
correlations in the CCF.

##  Appendix C Zero-padding

For a sequence of $N$ numbers $\\{x_{n}\\}$ in the time domain and
$\\{X_{k}\\}$ in the Fourier domain, the DFT simply follows

$X_{k}=\sum_{n=0}^{N-1}x_{n}\exp({-\frac{2\pi i}{N}kn}).$ (C1)

Zero-padding adds zeros to the end of the signal. Let $N$ be the intrinsic
number of inputs of the original discrete signal and $N^{\prime}$ be the total
number of inputs for the DFT with $(N^{\prime}-N)$ zeros padded after the
original signal. Denote the new signal as $\\{x^{\prime}_{n}\\}$. According to
Eq. C1, the Fourier transform of the zero-padded signal is

$\begin{split}X^{\prime}_{k}&=\sum_{n=0}^{N^{\prime}-1}x^{\prime}_{n}\exp({-\frac{2\pi
i}{N^{\prime}}kn})\\\ &=\sum_{n=0}^{N-1}x_{n}\exp({-\frac{2\pi
i}{N^{\prime}}kn})+\sum_{n=N}^{N^{\prime}-1}0\exp({-\frac{2\pi
i}{N^{\prime}}kn})\\\ &=\sum_{n=0}^{N-1}x_{n}\exp({-\frac{2\pi
i}{N^{\prime}}kn})\end{split}$ (C2)

Replace $k/N$ in Eq. C1 by $\xi$ and and $k/N^{\prime}$ Eq. C2 by
$\xi^{\prime}$, we have

$X_{k}=\sum_{n=0}^{N-1}x_{n}\exp({-2\pi i\xi n})$ (C3)

and

$X^{\prime}_{k}=\sum_{n=0}^{N-1}x_{n}\exp({-2\pi i\xi^{\prime}n}).$ (C4)

For $N^{\prime}>N$, $\xi^{\prime}$ means a finer sampling than $\xi$ in the
frequency grid. Therefore, the resulting amplitudes and phases in the Fourier
transform space look smoother, but it does not extract more features of the
original signal.

## References

  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252, doi: 10.1109/TPAMI.2015.2448083
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373
  * Breiman et al. (1984) Breiman, L., Friedman, J. H., Olshen, R. A., & Stone, C. J. 1984, Classification and Regression Trees (Monterey, CA: Wadsworth and Brooks)
  * Chaplin et al. (2019) Chaplin, W. J., Davies, G. R., Ball, W. H., Cegla, H. M., & Watson, C. A., E.-m. w. 2019, Astronomical Journal (New York, N.Y. Online), 157, doi: 10.3847/1538-3881/AB0C01
  * Collier Cameron et al. (2019) Collier Cameron, A., Mortier, A., Phillips, D., et al. 2019, MNRAS, 487, 1082, doi: 10.1093/mnras/stz1215
  * Collier Cameron et al. (2021) Collier Cameron, A., Ford, E. B., Shahaf, S., et al. 2021, MNRAS, 505, 1699, doi: 10.1093/mnras/stab1323
  * Cooley & Tukey (1965) Cooley, J., & Tukey, J. 1965, Mathematics of Computation, 19, 297
  * D’agostino et al. (1990) D’agostino, R. B., Belanger, A., & Jr., R. B. D. 1990, The American Statistician, 44, 316, doi: 10.1080/00031305.1990.10475751
  * de Beurs et al. (2020) de Beurs, Z. L., Vanderburg, A., Shallue, C. J., et al. 2020, Identifying Exoplanets with Deep Learning. IV. Removing Stellar Activity Signals from Radial Velocity Measurements Using Neural Networks.  https://arxiv.org/abs/2011.00003
  * Delchambre (2014) Delchambre, L. 2014, Monthly Notices of the Royal Astronomical Society, 446, 3545, doi: 10.1093/mnras/stu2219
  * Dumusque et al. (2014) Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132, doi: 10.1088/0004-637x/796/2/132
  * Dumusque et al. (2015) Dumusque, X., Glenday, A., Phillips, D. F., et al. 2015, The Astrophysical Journal, doi: 10.1088/2041-8205/814/2/l21
  * Dumusque et al. (2021) Dumusque, X., Cretignier, M., Sosnowska, D., et al. 2021, A&A, 648, A103, doi: 10.1051/0004-6361/202039350
  * Faria et al. (2022) Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115, doi: 10.1051/0004-6361/202142337
  * Fischer et al. (2016) Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, PASP, 128, 066001, doi: 10.1088/1538-3873/128/964/066001
  * Genton (2001) Genton, M. G. 2001, Journal of Machine Learning Research, 2, 299.  http://www.jmlr.org/papers/v2/genton01a.html
  * Gilbertson et al. (2020a) Gilbertson, C., Ford, E. B., & Dumusque, X. 2020a, Research Notes of the AAS, 4, 59, doi: 10.3847/2515-5172/ab8d44
  * Gilbertson et al. (2020b) Gilbertson, C., Ford, E. B., Jones, D. E., & Stenning, D. C. 2020b, The Astrophysical Journal, 905, 155, doi: 10.3847/1538-4357/abc627
  * Hastie (2015) Hastie, T. 2015, Statistical learning with sparsity : the lasso and generalizations, Chapman & Hall/CRC monographs on statistics & applied probability ; 143 (Boca Raton, FL: CRC Press)
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., & Friedman, J. 2009, The elements of statistical learning: data mining, inference and prediction, 2nd edn. (Springer).  http://www-stat.stanford.edu/~tibs/ElemStatLearn/
  * Heinisch (1962) Heinisch, O. 1962, Biometrische Zeitschrift, 4, 207
  * Herbst & Miller (1989) Herbst, W., & Miller, J. R. 1989, AJ, 97, 891, doi: 10.1086/115035
  * Jones et al. (2020) Jones, D. E., Stenning, D. C., Ford, E. B., et al. 2020, Improving Exoplanet Detection Power: Multivariate Gaussian Process Models for Stellar Activity.  https://arxiv.org/abs/1711.01318
  * Langellier et al. (2021) Langellier, N., Milbourne, T. W., Phillips, D. F., et al. 2021, AJ, 161, 287, doi: 10.3847/1538-3881/abf1e0
  * Matloch, L. et al. (2009) Matloch, L., Cameron, R., Schmitt, D., & Schüssler, M. 2009, A&A, 504, 1041, doi: 10.1051/0004-6361/200811200
  * Minasny & McBratney (2005) Minasny, B., & McBratney, A. B. 2005, Geoderma, 128, 192, doi: https://doi.org/10.1016/j.geoderma.2005.04.003
  * Pearson et al. (2018) Pearson, K. A., Palafox, L., & Griffith, C. A. 2018, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/stx2761
  * Pepe et al. (2002) Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632, doi: 10.1051/0004-6361:20020433
  * Pepe et al. (2014) Pepe, F., Molaro, P., Cristiani, S., et al. 2014, Astronomische Nachrichten, 335, 8, doi: https://doi.org/10.1002/asna.201312004
  * Petersburg et al. (2020) Petersburg, R. R., Ong, J. M. J., Zhao, L. L., et al. 2020, AJ, 159, 187, doi: 10.3847/1538-3881/ab7e31
  * Phillips et al. (2016) Phillips, D. F., Glenday, A., Dumusque, X., et al. 2016, Proceedings of SPIE, doi: 10.1117/12.2232452
  * Press et al. (2007) Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University Press).  https://books.google.com/books?id=1aAOdzK3FegC
  * Queloz et al. (2001) Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279, doi: 10.1051/0004-6361:20011308
  * Queloz et al. (2009) Queloz, D., Bouchy, F., Moutou, C., et al. 2009, A&A, 506, 303, doi: 10.1051/0004-6361/200913096
  * Rajpaul et al. (2015) Rajpaul, V. M., Aigrain, S., Aigrain, S., et al. 2015, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/stv1428
  * Raju et al. (1997) Raju, N. S., Bilgic, R., Edwards, J. E., & Fleer, P. F. 1997, Applied Psychological Measurement, 21, 291, doi: 10.1177/01466216970214001
  * Schwab et al. (2019) Schwab, C., Bender, C., Blake, C., et al. 2019, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 408.03
  * Shapiro & Wilk (1965) Shapiro, S. S., & Wilk, M. B. 1965, Biometrika, 52, 591, doi: 10.1093/biomet/52.3-4.591
  * Strassmeier, K. G. et al. (2018) Strassmeier, K. G., Ilyin, I., & Steffen, M. 2018, A&A, 612, A44, doi: 10.1051/0004-6361/201731631
  * Tibshirani (1996) Tibshirani, R. 1996, Journal of the Royal Statistical Society. Series B (Methodological), 58, 267.  http://www.jstor.org/stable/2346178
  * Vázquez Ramió, H. et al. (2005) Vázquez Ramió, H., Régulo, C., & Roca Cortés, T. 2005, A&A, 443, L11, doi: 10.1051/0004-6361:200500191
  * Zhao (2019) Zhao, J. 2019, PhD thesis, University of New South Wales, doi: 10.26190/unsworks/21571
  * Zhao (2022) Zhao, J. 2022, jinglinzhao/FIESTA-II: v2.1, v2.1, Zenodo, doi: 10.5281/zenodo.6554094
  * Zhao & Tinney (2020) Zhao, J., & Tinney, C. G. 2020, MNRAS, 491, 4131, doi: 10.1093/mnras/stz3254
  * Zhao et al. (2020) Zhao, L., Fischer, D. A., Ford, E. B., et al. 2020, Research Notes of the American Astronomical Society, 4, 156, doi: 10.3847/2515-5172/abb8d0
  * Zhao et al. (2022) Zhao, L. L., Fischer, D. A., Ford, E. B., et al. 2022, AJ, 163, 171, doi: 10.3847/1538-3881/ac5176