A statistical primer on exoplanet detection methods

Naomi Giertych, Jonathan P Williams, and Pragya Haravu
Department of Statistics, North Carolina State University

###### Abstract

Historically, a lack of cross-disciplinary communication has led to the
development of statistical methods for detecting exoplanets by astronomers,
independent of the contemporary statistical literature. The aim of our paper
is to investigate the properties of such methods. Many of these methods (both
transit- and radial velocity-based) have not been discussed by statisticians
despite their use in thousands of astronomical papers. Transit methods aim to
detect a planet by determining whether observations of a star contain a
periodic component. These methods tend to be overly rudimentary for starlight
data and lack robustness to model misspecification. Conversely, radial
velocity methods aim to detect planets by estimating the Doppler shift induced
by an orbiting companion on the spectrum of a star. Many such methods are
unable to detect Doppler shifts on the order of magnitude consistent with
Earth-sized planets around Sun-like stars. Modern radial velocity approaches
attempt to address this deficiency by adapting tools from contemporary
statistical research in functional data analysis, but more work is needed to
develop the statistical theory supporting the use of these models, to expand
these models for multiplanet systems, and to develop methods for detecting
ever smaller Doppler shifts in the presence of stellar activity.

Keywords: astrostatistics, Doppler shift, Gaussian process, radial velocity,
transit methods

##  1 Introduction

The discovery of an exoplanet orbiting the star 51 Pegasi (Mayor & Queloz,
1995) launched the exoplanet detection and characterization subfield of
astronomy. First and second generation spectrographs, however, could only
detect planets with larger radial velocities (i.e., planets with larger
planet-to-parent-star mass ratios) and not planets inducing a sub-1 m/s radial
velocity (RV). With third generation spectrographs, however, the limiting
factor for detecting Earth-like planets around Sun-like stars is not
instrument precision but rather deficiencies in statistical techniques for
decoupling signal from noise, even with potentially massive amounts of data.

The aim of our paper is to investigate and demonstrate the theoretical and
empirical properties of current and classical statistical methods used to
detect exoplanets. To the best of our knowledge, no such review has been
published from a statistician’s perspective, as many of the classical methods
were developed independently by astronomers. In particular, it is pointed out
in Feigelson et al. (2021), with reference to the classical phase dispersion
minimization (PDM, Stellingwerf, 1978) periodogram, that:

> “None of these statistical issues [relating to the PDM] have been examined
> by mathematical statisticians despite the use of the PDM in $\sim$1,500
> astronomical papers over four decades.”

Our paper contributes to filling this void as well as other voids in the
statistical review of exoplanet detection methods. The goal is to foster
further communications between the statistical and astronomical research
communities and to attract the attention of more statisticians to the young
and growing field of astrostatistics.

In the astronomical literature, statistical methods for exoplanet detection,
both classical and modern, can broadly be categorized as either transit-based
methods or RV-based methods. Many classical transit-based methods are
characterized as period detection techniques such as the phase-folding
strategies proposed by Stellingwerf (1978), Schwarzenberg-Czerny (1989),
Dworetsky (1983), and Lafler & Kinman (1965) or the Lomb-Scargle periodogram
proposed in Lomb (1976); Scargle (1982). For the RV approach, classical
methods include the cross-correlation function (CCF) introduced in Baranne et
al. (1996). More recent methods are based on Gaussian process (GP) regression
(e.g., Jones et al., 2017) or employ clever formulations of functional data
models inspired by extremely small scale of Doppler shifts (e.g., Holzer et
al., 2020). The functional data approach is particularly interesting because
the field of functional data analysis is still a rapidly developing area of
research in statistics communities. Further, merging the statistical field of
functional data analysis with exoplanet detection applications is critical for
the next generation of breakthroughs in astrostatistics methodologies
necessary with the 500 petabytes of data expected from the 10-year Legacy
Survey of Space and Time expected to begin in 2023.

In Section 2, we provide an in-depth examination of the classical transit-
based methods. The mathematical details are presented for how each method is
constructed, followed by discussions of the correctness of the ensuing
statistical inferences that were prescribed by the astronomy community for the
use of each method. In each case, we make comparisons to the appropriate
perspective from the statistical literature and provide extensive simulation
and empirical analyses to support our assessments. In general, we find the
classical transit-based methods to be overly rudimentary in their application
on the highly complex data they were designed for and to lack robustness to
model-specification assumptions. These are critical concerns because
exoplanetary data are well-known for their high degree of heteroscedasticity
and measurement errors. Next, in Section 3 we review stellar activity and the
basic astronomy of RV and Doppler shifts and then analyze classical and modern
RV-based approaches to exoplanet detection. The modern RV approaches tie into
contemporary statistical research in functional data analysis and naturally
lead to promising directions for future research.

##  2 Transit methods

Transit-based methods for exoplanet detection are developed on the assumption
that the observed amount of light a star emits (i.e., its flux) as a planet
orbits around it is a periodic function. With this assumption, the transit
method strategy for identifying exoplanets is to determine the statistical
significance of a nonzero period associated with star flux measurements over
time. Accordingly, the statistical model for an observed flux measurement (for
a given star) at time $t\in\\{t_{1},\dots,t_{T}\\}$ is

$X_{t}=h(t)+\epsilon_{t},$ (1)

where $h(t)$ is a periodic function with true period $p_{0}$ (i.e.,
$h(t)=h(t+kp_{0})$, for any integer $k$), and
$\epsilon_{t_{1}},\dots,\epsilon_{t_{T}}\overset{\text{iid}}{\sim}N(0,\sigma^{2})$
for some $\sigma>0$. Classical statistical inference developed by astronomers,
using data arising from this model are predominantly either phase-based or
based on Fourier transforms. These methods are described in the following two
sections.

###  2.1 Phase-based methods

Figure 1: The leftmost panel displays synthetic flux values generated as
perturbations from a periodic function. The middle and rightmost panels
illustrate the result of phase-folding the flux values according to an
incorrect and the correct period, respectively.

Phase-based methods aim to determine the period of a periodic function based
on noisy observations, and they leverage the idea that observations can be
folded (or phase-folded) in time modulo the true period of the function. Given
a trial period, observed flux data can always be phase-folded, but if the
trial period is not close to the true period then the phase-folded values
likely will not align with the function $h$ in phase; see Figure 1. In the
classical literature on these methods, Stellingwerf (1978) and Schwarzenberg-
Czerny (1989) suggest binning the phase-folded data into groups with similar
phase-values and deriving statistics based on ratios of sums of squared
deviations to test various trial periods. Namely, for a trial period $p$, the
null and alternative hypotheses suggested by Stellingwerf (1978) are
formulated as

$\begin{split}H_{0}&:\mu_{1}=\cdots=\mu_{r}\hskip 18.06749pt(\text{i.e.,
either }p\neq p_{0}\text{ or }p_{0}=\infty)\\\ H_{a}&:\text{at least one
$\mu_{i}$ is different from the others}\hskip 18.06749pt(\text{i.e.,
}p=p_{0}\text{ or }p_{0}<\infty)\\\ \end{split}$ (2)

where $\mu_{1},\dots,\mu_{r}$ are defined as the bin means for the phase-
folded data modulo $p$. The bin means being equal (as in $H_{0}$) suggests no
periodicity in the observed flux values, and so statistical evidence
supporting the assertion that a planet is orbiting the star is established if
a test rejects $H_{0}$ in favor of $H_{a}$. Two such classical tests are
constructed, next.

Let $X_{i,j}$ be an observed flux value from a time series which has already
been binned in-phase based on some trial period $p$, where
$i\in\\{1,\dots,r\\}$ is the bin index, $r$ is some specified number of bins,
$j\in\\{1,\dots n_{i}\\}$ is the observation index for the $i$-th bin, and
$n_{i}$ is the number of observations in bin $i$. Let $n:=\sum_{i}n_{i}$ be
the total number of observed time instances, and take $S_{0}^{2}$,
$S_{1}^{2}$, and $S_{2}^{2}$ to denote the total sample variation, the across-
bins mean variation, and the within-bin variation, respectively. Then,

$\begin{split}(n-1)S_{0}^{2}&=\sum_{i=1}^{r}\sum_{j=1}^{n_{i}}(X_{i,j}-\bar{X})^{2}\\\
&=\sum_{i=1}^{r}n_{i}(\bar{X}_{i}-\bar{X})^{2}+\sum_{i=1}^{r}\sum_{j=1}^{n_{i}}(X_{i,j}-\bar{X_{i}})^{2}=(r-1)S_{1}^{2}+(n-r)S_{2}^{2}.\end{split}$
(3)

The PDM statistic proposed by Stellingwerf (1978) and the analysis of variance
(AOV) statistic proposed by Schwarzenberg-Czerny (1989) are defined,
respectively, as

$\displaystyle T_{\text{PDM}}$ $\displaystyle:=S_{2}^{2}/S_{0}^{2}$ (4)
$\displaystyle T_{\text{AOV}}$ $\displaystyle:=S_{1}^{2}/S_{2}^{2},$ (5)

and the trial period $p$ is tested by determining the tail probability of the
sampling distributions of these statistics under the null hypothesis $H_{0}$.

Under $H_{0}$ and assuming independently and identically Gaussian distributed
noise in the observed flux values, Stellingwerf (1978) incorrectly argued that
the PDM statistic followed some parameterization of the F distribution.
Unfortunately, there has been much confusion about these facts in the
astronomy literature, largely due to a lack of commentary from the statistics
community. The inappropriateness of the F distribution for the PDM statistic
is discussed in Heck et al. (1985) and referenced in Linnell Nemec & Nemec
(1985) with incorrect mathematical justification for why it does not follow an
F distribution. These authors incorrectly claim that the PDM statistic does
not follow an F distribution because the numerator, $S_{2}^{2}$, does not
follow a chi-squared distribution. They are mislead in remarking that the
matrix associated with its quadratic form is not idempotent. Subsequently,
however, Schwarzenberg-Czerny (1989) correctly concluded that the PDM
statistic does not follow a F distribution since the numerator and denominator
are not independent (as is clear from equation (3)). In a followup paper,
Schwarzenberg-Czerny (1997) correctly argues that a scaled version of the PDM
statistic, specifically $\frac{n-r}{n-1}T_{\text{PDM}}$, follows the beta
distribution with shape parameters $\frac{n-r}{2}$ and $\frac{r-1}{2}$,
assuming $H_{0}$ is true. Additionally, Schwarzenberg-Czerny (1989) correctly
states that under $H_{0}$ the AOV statistic follows the F distribution with
degrees of freedom $r-1$ and $n-r$.

Arguing that $T_{\text{AOV}}\sim\text{F}_{r-1,n-r}$ is a direct application of
Cochran’s theorem after noting that equation (3) can be expressed in terms of
the quadratic forms

$X^{\prime}(I_{n}-n^{-1}J_{n})X=X^{\prime}(P_{V}-n^{-1}J_{n})X+X^{\prime}(I_{n}-P_{V})X,$
(6)

where $X:=(X_{11},\dots,X_{1n_{1}},\dots,X_{r1},\dots,X_{rn_{r}})^{\prime}$,
$I_{n}$ is the $n\times n$ identity matrix, $J_{n}$ is an $n\times n$ matrix
of ones, and $P_{V}:=V(V^{\prime}V)^{-1}V^{\prime}$ is the projection matrix
onto the column space of $V:=\text{diag}\\{1_{n_{1}},\dots,1_{n_{r}}\\}$
corresponding to the design matrix in the linear regression $X=V\mu+\epsilon$
with $\mu=(\mu_{1},\dots,\mu_{r})^{\prime}$ and
$\epsilon:=(\epsilon_{11},\dots,\epsilon_{1n_{1}},\dots,\epsilon_{r1},\dots,\epsilon_{rn_{r}})^{\prime}$.
Note that $V$ is simply the design matrix associated with the phase-folding of
the original observed flux values modulo the trial period $p$. Cochran’s
theorem applies to establish that $S_{1}^{2}/\sigma^{2}\sim\chi_{r-1}^{2}$
independently of $S_{2}^{2}/\sigma^{2}\sim\chi_{n-r}^{2}$ because from
equation (6),

$X^{\prime}X=X^{\prime}\underbrace{n^{-1}J_{n}}_{=:A_{1}}X+X^{\prime}\underbrace{(P_{V}-n^{-1}J_{n})}_{=:A_{2}}X+X^{\prime}\underbrace{(I_{n}-P_{V})}_{=:A_{3}}X,$

with $A_{1},A_{2}$, and $A_{3}$ all idempotent, $A_{1}+A_{2}+A_{3}=I_{n}$, and
rank$(A_{1})=1$, rank$(A_{2})=r-1$, and rank$(A_{3})=n-r$.

Arguing that
$T_{\text{PDM}}\sim\text{beta}\big{(}\frac{n-r}{2},\frac{r-1}{2}\big{)}$
follows by a series of arithmetic steps to obtain

$\begin{split}F_{\frac{n-r}{n-1}T_{\text{PDM}}}(u)&=P\Big{(}\frac{n-r}{n-1}T_{\text{PDM}}\leq
u\Big{)}=1-F_{T_{\text{AOV}}}\Big{(}\frac{n-r}{r-1}(u^{-1}-1)\Big{)}.\\\
\end{split}$

Then, since $T_{\text{AOV}}\sim\text{F}_{r-1,n-r}$, taking derivatives on both
sides yields, for $u\in[0,1]$,

$f_{\frac{n-r}{n-1}T_{\text{PDM}}}(u)=\frac{u^{\frac{n-r}{2}-1}(1-u)^{\frac{r-1}{2}-1}}{B(\frac{n-r}{2},\frac{r-1}{2})};$

the density function of the
$\text{beta}\big{(}\frac{n-r}{2},\frac{r-1}{2}\big{)}$.

Observe that the $\sigma^{2}$ terms cancel when constructing $T_{\text{AOV}}$
and $T_{\text{PDM}}$. A point of caution is that this cancellation will not
occur in the case that the errors are heteroscedastic, invalidating the
appropriateness of relying on critical values from the F or beta distribution,
respectively. In the heteroscedastic setting, the practitioner would need to
know the variances to scale each observation to correct the misspecification.

If the alternative hypothesis $H_{a}$ is true (i.e., if there is a planet
orbiting the star) the distributions of $T_{\text{AOV}}$ and $T_{\text{PDM}}$
depend on whether the trial period $p$ is the true period $p_{0}$ or not. In
the case that $p=p_{0}$, it is established from well-known statistical theory
that $T_{\text{AOV}}$ follows the non-central F distribution
$\text{F}_{r-1,n-r}\big{(}\frac{1}{2\sigma^{2}}\mu^{\prime}A_{2}\mu\big{)}$
and $T_{\text{PDM}}$ follows the Type II non-central beta distribution
$\text{beta}_{\frac{n-r}{2},\frac{r-1}{2}}\big{(}\frac{1}{2\sigma^{2}}\mu^{\prime}A_{2}\mu\big{)}$.
In the case that $p\neq p_{0}$, the flux values have been folded modulo the
wrong period, and so the true data generation equation has the form
$X=W\mu+\epsilon$ for some design matrix $W\neq V$. Accordingly,
$T_{\text{AOV}}$ and $T_{\text{PDM}}$ follow doubly non-central F and doubly
non-central beta distributions, respectively. While the distributions that
arise under $H_{a}$ could theoretically be used to determine the power of the
tests, in practice the values of $\mu$ and $\sigma$ are not known.

Additionally, in the context of exoplanet detection, since the true period of
the exoplanet is unknown, many trial periods might be tested. These tests
amount to constructing a confidence interval around the true period where the
period associated with the smallest PDM statistic (or largest AOV statistic)
might be used to as an estimate of the true period; e.g., see Figure 8.

Lastly, it is worth commenting on a more rudimentary but also popular phase-
based method, the string-length method, which pre-dates the sums of squares
based analyses. Phase-folding observed flux data values based on a trial
period $p$ is defined by the phase-scale and flux-value pair
$(\phi_{t},X_{t})$ at time $t$, where $\phi_{t}:=(t\mod p)/p$ and $X_{t}$ is
generated as in equation (1). It is suggested in Burke Jr et al. (1970) to
measure the fit of the trial period $p$ based on the Euclidean distance
between neighboring time-flux pairs after phase-folding; since the phase-
folded values are periodic, the distance between the first and last points is
also included with an adjustment for the difference in phase values.
Specifically, the string-length is defined

$L:=[(X_{t_{1}}-X_{t_{T}})^{2}+(\phi_{t_{1}}-\phi_{t_{T}}+1)^{2}]^{1/2}+\sum_{t=t_{1}}^{t_{T}}[(X_{t}-X_{t-1})^{2}+(\phi_{t}-\phi_{t-1})^{2}]^{1/2}.$

The period $p$ that minimizes this distance is taken as the estimated period
for the observed data values. Dworetsky (1983) examines the statistical
properties of this distance and obtains estimates for the population mean
length and standard deviation; however, he only considers circular orbits and
Gaussian distributed errors for these estimates. Additionally, while the
string-length method avoids the arbitrary binning of observations used in the
PDM and AOV statistics, Feigelson et al. (2021) notes that it is less
sensitive to detection of the true period in practice.

###  2.2 Lomb-Scargle periodogram

The phase-based methods discussed in the previous section are designed to work
for time domain values of periodic functions. Methods based on Fourier
transforms have also been considered, however, for analysis in the frequency
domain. The intuition is that the true period of a periodic function will be
exhibited as a spike in its frequency spectrum.

Assuming the time domain flux values are generated as in equation (1), the
Fourier transform is defined as

$\hat{X}(f):=\int_{-\infty}^{\infty}X_{t}e^{-2\pi ift}dt,$

for frequency $f\in\mathbb{R}$. The square of the Fourier transform is
referred to as the power spectral density and is a measure of how much each
frequency characterizes a periodic function (Priestley, 1981; VanderPlas,
2018). The Lomb-Scargle periodogram aims to estimate the power spectral
density, particularly for the case of exoplanet data with unevenly sampled
observations in time (i.e., when classical time series methods do not apply).
The Lomb-Scargle periodogram for trial frequency $f$ is

$P(f):=\frac{1}{2}\bigg{(}\frac{\\{\sum_{t=t_{1}}^{t_{T}}X_{t}\cos[2\pi
f(t-\tau)]\\}^{2}}{\sum_{t=t_{1}}^{t_{T}}\cos^{2}[2\pi
f(t-\tau)]}+\frac{\\{\sum_{t=t_{1}}^{t_{T}}X_{t}\sin[2\pi
f(t-\tau)]\\}^{2}}{\sum_{t=t_{1}}^{t_{T}}\sin^{2}[2\pi f(t-\tau)]}\bigg{)},$

where $\tau:=\frac{1}{4\pi
f}\tan^{-1}\bigg{[}\frac{\sum_{t=t_{1}}^{t_{T}}\sin(4\pi
ft)}{\sum_{t=t_{1}}^{t_{T}}\cos(4\pi ft)}\bigg{]}$ (VanderPlas, 2018). When
the data are generated as Gaussian noise, $2P(f)/\sigma^{2}\sim\chi_{2}^{2}$
for each $f$, established by Cochran’s theorem and the fact that
$\sum_{t}\sin[2\pi f(t-\tau)]\cos[2\pi f(t-\tau)]=0$, for any $f$.

Figure 2: Comparison of synthetic Lomb-Scargle Periodogram false alarm
probabilities versus the theoretical bound provided by Baluev (2008).

One of the main difficulties of period detection using the Lomb-Scargle
periodogram is determining the significance of the largest peak. This
difficulty is due to the periodogram value at one frequency being correlated
with the value at other frequencies through a convolution of the true signal
and the observation times, where a convolution is defined as
$f*g(t):=\int_{\infty}^{\infty}f(\tau)g(t-\tau)d\tau$ (VanderPlas, 2018). In
the case of Gaussian errors as in equation (1), Baluev (2008) adapted theory
(based on work from Davies (1977), Davies (1987), and Davies (2002)) for
approximating the distribution of the maximum of chi-squared random variables
over frequencies $f$, for the Lomb-Scargle periodogram peak. While this
approximation gives an analytic result that can be easily implemented, it is
an upper bound for periodograms and can greatly overestimate the type 1 error
rate when aliasing occurs (see Section 2.3 and VanderPlas, 2018).

In Figure 2, we examine the type 1 error rate of the Lomb-Scargle periodogram
compared to the upper bound derived in Baluev (2008). We follow VanderPlas
(2018) and generate 100 unevenly sampled observations $m_{1},\dots,m_{100}$
from the set $\\{0,1,\dots,5\cdot 24\\}$ (i.e., 100 observations times over 5
days) and rescale the observation times as $t_{i}=0.01\cdot m_{i}$ for
$i\in\\{1,\dots,100\\}$. Following VanderPlas (2018), we generate a grid of
frequency values from $1/(24\cdot 25)$ to $24\cdot 10$. Finally, for each of
1,000 iterations, we randomly generate the flux values from a standard
Gaussian distribution (i.e., under the null hypothesis) and calculate the
Lomb-Scargle power at each point in the frequency grid. Finally, the empirical
false alarm probabilities (i.e., p-values) are calculated (using the
`LombScargle` function from the Python module Astropy (Robitaille et al.,
2013)) for the range of maximum periodogram values observed over the 1,000
iterations. These empirical false alarm probabilities are plotted in Figure 2,
along with the bound from Baluev (2008).

###  2.3 Overarching limitations of transit methods

A primary difficulty encountered with transit methods is the limitation in the
range of periods that are identifiable for a periodic function, given the
denseness of available observation times. A periodic function with a given
frequency (the inverse of the period) must be sampled at a rate that is at
least twice the frequency in order to correctly reconstruct the signal; this
rate is known as the Nyquist rate (Leis, 2011). If the signal is sampled at a
rate lower than the Nyquist rate, then aliasing of the signal occurs,
resulting in observing a signal with a frequency shorter (period longer) than
the true frequency (Priestley, 1981). We illustrate this phenomenon in Figure
3 for a cosine curve observed at two sampling rates.

Figure 3: A sine curve, such as the solid green line, can be correctly
reconstructed by sampling observations at equal time intervals at least as
short as the Nyquist rate; these samples correspond to the red squares.
However, the observation times corresponding to the blue dots fit equally well
to the dashed purple line as they do to the original sine curve. This
phenomenon is referred to as aliasing.

Given a set of evenly-spaced observations in time, the Nyquist frequency/
Nyquist limit is the highest frequency (shortest period) in the data that we
can hope to detect without aliasing (Leis, 2011; VanderPlas, 2018). Moreover,
while these are well-defined notions in the evenly-spaced sampling scenario,
the Nyquist frequency for unevenly-spaced observations is only defined when
all observation pairs have ratios that are rational (see Eyer & Bartholdi,
1999; VanderPlas, 2018). Though, even if the unevenly-spaced Nyquist frequency
is defined, it can be prohibitively large for practical use, and a reasonable
limit must be imposed based on previous knowledge (VanderPlas, 2018).

Other difficulties with transit methods arise from a disparity in the
assumptions made for theory versus those that are reasonable in practice.
While the computation of the discussed transit method statistics are
straightforward, the choice of a grid of trial period or frequency values to
consider is nontrivial, and a poor choice can result in entirely missing a
true period or frequency of an exoplanet orbit (VanderPlas, 2018). Even with a
well-choosen frequency grid, aliasing may still be an issue, such that the
periodograms of the Lomb-Scargle, PDM, and AOV statistics cannot distinguish
between the true period and its aliases. Additionally, the Lomb-Scargle
periodogram and phase-based methods are only used to test for one strong
signal (i.e., the orbit of one exoplanet in the data), but many known systems
have multiple planets. Attempts have been made to remove a detected signal
from data and reapply the Lomb-Scargle periodogram, but since signals are
generally not orthogonal, these techniques introduce problems for detecting
additional periods (Khan et al., 2017). Furthermore, all of these methods rely
on independent error assumptions, even though it is known that starlight often
has correlated noise (Khan et al., 2017). In fact, Feigelson et al. (2021)
argues “the most insidious problem is that many stars show autocorrelated but
aperiodic variations that interact with the irregular [observation times] to
produce false peaks in the periodogram for short duration data sets.” These
false peaks can lead to claims of a planet being discovered when, in fact, one
has not.

###  2.4 Empirical study of transit methods

Time | Phase-Scale | Bin | $\mu_{i}$ | Flux ($X_{i,j}$)
---|---|---|---|---
443.510962 | 0.116228 | 1 | -0.002118 | -0.000797
443.531395 | 0.117800 | 1 | -0.002118 | -0.00256
⋮ | ⋮ | ⋮ | ⋮ | ⋮
538.121136 | 0.393934 | 3 | -0.000308 | -0.000257
538.162004 | 0.397077 | 3 | -0.000308 | -0.001744
Table 1: Organization of the data generated from the Kepler-1520b time stamp.

We use publicly available observations of the exoplanet Kepler-1520b from
about 443 to 538 days after the start of the Kepler mission; these data are
available in the Python module Lightkurve (Lightkurve Collaboration et al.,
2018). During this time, Kepler-1520b is observed 4485 times with an irregular
observation cadence of about half-an-hour on average. Each observation has a
corresponding flux value that has been corrected for systematic artifacts. We
use these observations to construct data with realistic observation cadences
but known true flux values to study the properties of the PDM, AOV, and Lomb-
Scargle statistics. Additionally, we explore the effect of irregular
observations and correlated noise on the sensitivity for detecting an
exoplanet.

####  2.4.1 Sampling distributions of the AOV and PDM statistics

Figure 4: The first row of panels display the synthetic phase-folded flux
values generated as independent Gaussian realizations from 10 bins with
distinct means, at a similar scale as the flux values observed for the
Kepler-1520b data. The true bin means (calculated after phase-folding) are
indicated by the black lines. The second and third rows display histograms of
the AOV and PDM statistics, respectively. The solid black line is the density
of the central F distribution in the second row and the density of the central
beta distribution in the third row. The true period is $p_{0}=13$.

In this section, we study the empirical sampling distributions of the AOV and
PDM statistics in three distinct scenarios. In all three scenarios, we
consider synthetic periodic data that have been phase-folded modulo the
correct period, as well as incorrect periods. Since both the AOV and PDM
statistics, along with their respective sampling distributions (under the
hypotheses in equation (2)), are derived from classical sums of squares
decompositions of iid Gaussian data, the first scenario we consider is based
on synthetic independent Gaussian data generated with different means across
10 bins. The bins are constructed by phase-folding the real Kepler-1520b data
time stamps, as in Table 1, according to an arbitrary trial period (we choose
a period of 13 days) resulting in a phase-scale value in $(0,1)$, and then
defining the $i$-th bin as associated to phase-scale values in
$(\frac{i-1}{10},\frac{i}{10})$ for $i\in\\{1,2,\dots,10\\}$.

Figure 5: Same as Figure 4 except the data are generated using true values
from a fitted sine curve to the observed Kepler data and homoscedastic
Gaussian noise with a standard deviation level similar to the observed Kepler
data. The true period is $p_{0}=5.976$.

Next, we randomly assign the bin means as
$\mu_{i}\sim\text{uniform}(-.0025,.0025)$ and generate flux values
$X_{i,1},\dots,X_{i,n_{i}}\overset{\text{iid}}{\sim}\text{N}(\mu_{i},s^{2})$
for $i\in\\{1,\dots,10\\}$, where $s=.0015$ is the sample standard deviation
(to 4 decimal places) of the actual Kepler-1520b flux values. For a single set
of values
$\mu_{1},\dots,\mu_{10}\overset{\text{iid}}{\sim}\text{uniform}(-.0025,.0025)$,
we carry out a simulation study by generating 300 synthetic data sets as
described, each time computing the AOV and PDM statistics as in equations (4)
and (5) under a trial period, and comparing the resulting sampling
distributions to their respective theorized distributions under $H_{0}$. We
consider the empirical sampling distributions for trial periods that do and do
not match the true period used to generate the data. See Figure 4.

The major conclusion from Figure 4 is that phase-folding flux values using an
incorrect trial period will not necessarily result in bins with a common mean,
implying that the null hypothesis $H_{0}:\mu_{1}=\dots=\mu_{r}$ may be
rejected even if the trial period $p$ is not equal to the true period $p_{0}$.
This phenomenon is illustrated for the sampling distributions that arise after
phase-folding modulo the (incorrect) trial periods of both 9 and 12 but was
also observed in a range of other trial periods that we considered. Thus,
$H_{0}$ may be rejected for a reason other than due to chance (controlled by
the type 1 error probability) or the detection of a period in the data,
potentially leading to the false discovery of an exoplanet.

Figure 6: Same as Figure 4 except the data are generated using true values
from a fitted sine curve to the observed Kepler data and heteroscedastic
Gaussian noise. The true period is $p_{0}=5.976$.

Note that we repeated this analysis two additional times with the sample
standard deviations of the synthetic flux values set as $\frac{2}{3}\cdot s$
and $\frac{4}{3}\cdot s$, respectively, but the results were largely
unchanged. We also repeated the analyses using bin widths of .05
(corresponding to 20 bins) and .2 (corresponding to 5 bins), but the analyses
were, again, largely unchanged.

Next, for our second of three scenarios for studying the empirical sampling
distributions of the AOV and PDM statistics, we repeat the analysis of the
first scenario, but for synthetic sinusoidal data, considering both
homoscedastic and heteroscedastic errors. To make our analysis most
appropriate, we first fit a sine curve to the Kepler-1520b data using the
least squares optimize function `least_squares` from the SciPy module in the
Python programming language (Virtanen et al., 2020) and then generate
synthetic flux values as Gaussian perturbations on the fitted sine curve at
the observed time stamps for the Kepler-1520b data set. Homoscedastic
perturbations are drawn independently from the N$(0,s^{2})$ distribution,
where $s=.0015$ as in the first simulation scenario, and the results are
presented in Figure 5. Heteroscedastic perturbations are drawn independently
for the phase-folded values (using the true period) in bin $i$, for
$i\in\\{1,\dots,10\\}$, from the N$(0,\rho_{i}^{2})$ distribution with
$\rho_{i}\sim\text{lognormal}(\log(s)-.5^{2}/2,.5^{2})$ distribution, so that
$E(\rho_{i})=s$. The results corresponding to the heteroscedastic
perturbations are presented in Figure 6.

Figure 7: Same as Figure 4 except the data are a sample of 300 observed
Kepler-1520b points taken after 80% of the original data has been randomly
removed. The true period is $p_{0}=0.654$.

In our final of three scenarios for studying the empirical sampling
distributions of the AOV and PDM statistics, we repeat the analysis of the
first two scenarios, but this time subsampling from the real flux values
observed in the Kepler-1520b data set. We calculate the AOV and PDM statistics
after phase-folding 300 randomly subsampled time-flux pairs according to the
three trial periods .640, .651, and .654 days. Note that it is widely agreed
upon in the astronomy community that the true period of Kepler-1520b is
approximately .654 days (exoplanets.org). To study the effect of various
missing data schemes that are common to exoplanet data sets, we consider five
different settings: (i) almost evenly-spaced observations (original data),
(ii) randomly omit 80 percent of the data prior to subsampling, (iii) remove 8
consecutive hours from each recorded day prior to subsampling, (iv) prior to
subsampling remove 16 consecutive hours from each recorded day, and (v) remove
2 consecutive days from each recorded week prior to subsampling.

Since the results did not differ much between scenarios, we only show the
subsampled and phase-folded flux values and the corresponding sampling
distributions of the AOV and PDM statistics under scenario (ii) in Figure 7.
Figure 7 again highlights that the empirical distributions of $T_{\text{AOV}}$
and $T_{\text{PDM}}$ may not follow central F and beta distributions,
respectively, under an incorrect trial period (i.e., when $H_{0}$ is should to
be true).

####  2.4.2 AOV, PDM, and Lomb-Scargle Periodograms

Since astronomers need to compute the PDM, AOV, and Lomb-Scargle values for
multiple trial periods/frequencies, it is often helpful to construct a
periodogram for associating the value of these statistics against their trial
periods/frequencies. Periodograms, thus, help determine significant trial
periods/frequencies. In this section, we examine how these periodograms
visually change under five scenarios: (a) flux data generated from an evenly-
sampled sine curve and iid noise, (b) flux data generated from a randomly-
sampled sine curve and iid noise, (c) flux data generated from an evenly-
sampled sine curve during weeknights and iid noise, (d) flux data generated
from an evenly-sampled saw-tooth curve and iid noise, and (e) flux data
generated from an evenly-sampled sine curve and correlated noise. These
scenarios are simplifications of the data generation schemes astronomers
contend with. Scenario (a) is the baseline with which to compare all other
scenarios.

Figure 8: The columns of panels correspond to periodograms over a range of
trial periods for the AOV, PDM, and Lomb-Scargle power statistics,
respectively. The first three rows present the periodograms for iid
perturbations from a sine curve with equal observation times, random
observation times, and observation only during weeknights, respectively (i.e.,
scenarios (a), (b), and (c)). The fourth row presents the periodograms for iid
perturbations from a saw-tooth curve with equal observation times (i.e.,
scenario (d)), and the last row corresponds to iid perturbations from a sine
curve with equal observation times and correlated noise (i.e., scenario (e)).
The solid black line is the mean periodogram for 1,000 simulated data sets, at
each trial period, and the dotted blue lines are the mean periodogram plus and
minus the Monte Carlo standard deviation, at each trial period. The true
period is indicated by a vertical red line.

Scenario (a) is generated using a period of $2\pi$ and is sampled at a step
size of 0.0211 until 4,486 observations (similar to Kepler-1520b, above) have
been made; finally iid Gaussian noise with a standard deviation of 0.0016
(again, similar to Kepler-1520b) is added to the data. Scenario (b) is
generated by randomly removing 85 percent of the data, a level expected from
astronomical data based on a talk given at the 2021 Statistical Challenges in
Modern Astronomy VII conference. Scenario (d) is generated in the same way as
scenario (a) but with the sine curve replaced by a saw-tooth curve. Finally,
scenario (e) is generated by replacing the iid Gaussian noise from scenario
(a) with correlated noise using equation 11 as the covariance kernel. The
covariance kernel parameters are set such that the diagonal components are the
same as in the iid case, and observations within $\pm 5$ neighboring time
points are correlated.

We generate 1,000 synthetic data sets for each of the scenarios above. For
each synthetic data set, the AOV and PDM statistics are calculated for trial
periods ranging from $\pi$ to $3\pi$ with a step size of 0.00628, and the
Lomb-Scargle power is estimated using the default frequency grid in the Python
module Astropy (Robitaille et al., 2013). Figure 8 shows the mean observed
AOV, PDM, and Lomb-Scargle periodograms. Examining scenario (b), it appears
that all of the statistics are robust to missing completely at random
observations in time, while scenario (c) suggests that the AOV and PDM
statistics periodograms are corrupted from systematically missing observation
times. The results of panel (c) are particularly worrisome due to the shared
scheduling of telescope resources that astronomers must navigate when
requesting collection of data on a potentially undiscovered exoplanet. All
three statistics appear to be robust to non-sinusoidal signals in that they
all have a peak at the true period in scenario (d), though, the maximum of the
Lomb-Scargle periodogram also has peaks corresponding to aliases. Further
investigation is needed to separate the aliases from the true period, in
practice. None of the statistics, however, suggest evidence of robustness to
the correlated noise; such a finding is highly troublesome if, for example,
measurement errors are somehow correlated with flux values.

##  3 Radial velocity methods

Figure 9: Doppler shift of a star’s spectrum when the star has an orbiting
companion planet.

RV methods are an increasingly more popular approach than transit-based
methods for exoplanet detection. To introduce these RV-based methods we first
need to explain some necessary topics in astronomy. Light waves are measured
in fractions of centimeters called angstroms (Å) and recall that the flux of a
star is a measure of the amount of light it emits. The spectrum of a star is
the amount of flux the star emits at a collection of wavelengths ($\lambda$)
measured in angstroms. A lone star with no objects orbiting it has a constant
light spectrum over time. This constant light spectrum contains absorption
features where the amount of flux drops below the surrounding flux values (the
continuum) for a specific wavelength. A star with a companion planet, however,
orbits around a common center of mass (with the companion planet), called a
barycenter, located some distance away from the center of the star.

Figure 10: Spectrum of a star without a companion planet orbiting (solid,
black line), and Doppler shifts of the spectrum as the star orbits toward
(blue, dotted line) and away (red, dashed line) from Earth.

From our perspective on Earth, a star wobbles toward and away from us as it
orbits its barycenter. As the star moves toward (away from) Earth, the light
waves it emits are compressed (stretched), and the spectrum shifts toward the
blue (red) end. A shift in the spectrum is called a Doppler effect, and the
amount of shift is called the RV; a schematic of this effect is given in
Figure 9. A comparison of the true spectrum of the star with its blue- and
red-shifted spectra for a small RV is shown in Figure 10. RV-based methods for
exoplanet detection are essentially strategies for estimating and determining
the statistical significance of the RV from observed flux measurements of a
star.

Objects with negative RV are moving toward the observer, and objects with
positive RV are moving away from the observer. The Doppler effect at time $t$
is the wavelength rescaling that occurs due to the RV of an object with
respect to a fixed point, measured by

$\xi_{t}:=\frac{1+v_{r}(t)/c}{\sqrt{1-v^{2}/c^{2}}},$

where $v_{r}(t)$ is the RV, $c$ is the speed of light, and $v$ is the absolute
speed of the source (Holzer et al., 2020). If $f(\lambda)$ denotes the
normalized flux as a function of the wavelength $\lambda$, then
$f(\lambda\xi_{t})$ is the normalized flux of the Doppler shifted light. It is
argued in Holzer et al. (2020) that, assuming orbiting companions are moving
well below the speed of light and barycentric corrections have been made,
$\xi_{t}$ can reasonably be approximated as

$\xi_{t}=1+v_{r}(t)/c.$ (7)

The RV induced by a single planet system at time $t$ is given by

$v_{r}(t)=K(\zeta\cos(\omega)+\cos(\omega+\phi(t)))+\gamma,$ (8)

where $\zeta$ is the eccentricity of the planetary orbit, $\omega$ is the
argument of periapsis, $K$ is the amplitude, $\gamma$ is the velocity offset,
and $\phi(t)$ is the phase of the star in its elliptical orbit around the
barycenter (Jones et al., 2017). The argument of periapsis is the angle
between the black lines displayed in Figure 9.

###  3.1 Data collection and preprocessing

The observed spectrum of a star is measured with a spectrograph which
disperses light into its component wavelengths, much like a prism disperses
white light into the rainbow. The light is focused on a detector, called a
charged coupled device, that is composed of a grid of independent semi-
conductors, known as pixels (Baranne et al., 1996). While these span a range
of wavelengths, only the central value is recorded as the observed wavelength.
Specifically, the number of photons that hit a pixel during an observation
period is defined as the flux observed at the corresponding wavelength. In a
final step, the spectrograph collects the wavelength-flux pairs (i.e., the
observed spectrum) in partially overlapping groups based on wavelength values;
these groups are typically referred to as (wavelength) orders.

Researchers usually perform a series of preprocessing steps on raw flux values
before estimating the RV; we highlight the major steps, some of which affect
exoplanet detection techniques. First, the wavelength orders are usually
stitched together using an interpolation scheme (since the overlapping
wavelength-flux pairs are not exactly equal). Second, the spectrum is
corrected for instrument and constant astrophysical effects; the corrected
spectrum is what is often referred to as the normalized flux or normalized
spectrum. Third, the movement of the Earth around the Sun is accounted for in
what is referred to as the barycentric-correction. Lastly, for ground-based
telescopes, the star’s spectrum is contaminated by absorption features caused
by Earth’s atmosphere called tellurics. These tellurics are traditionally
removed from the spectrum, and we assume these corrections have been made
unless otherwise specified.

The current biggest challenge for the detection of small RVs is that stellar
activity can mimic a Doppler shift. Stellar magnetic fields can induce active
regions that, combined with the star’s own rotation, can lead to a spurious
signal with a period on the order of days to weeks (Rajpaul et al., 2015;
Jones et al., 2017). Additionally, long-term magnetic cycles can lead to RV
perturbations on the scale of a planet and with periods on the order of years
(Rajpaul et al., 2015). The sum of a RV and any corresponding stellar nuisance
signals is referred to as the corrupted RV. Fortunately, stellar activity
effects each absorption feature differently (e.g., by changing its shape),
whereas a Doppler shift effects the entire spectrum equally. Exploiting this
phenomenon has not yet been widely discussed in the literature.

While transit methods aim to identify a statistically significant period in
observed flux data, RV methods attempt to compare the spectrum of a star to
its Doppler shifted spectrum. Some methods aim to estimate the change in RV
between observed spectra (i.e., the RV that is needed to shift one observed
spectra to match another). Other methods endeavor to estimate the RV with
reference to the spectrum of the star without a planet (e.g., the black line
in Figure 10) by constructing a template mask (e.g. the circles in Figure 10).
The template mask is constructed either from knowledge of the type of star or
from observations of the star, and it represents the spectrum of the star
without a planet. A representative collection of RV-based approaches are
introduced and examined in the sections that follow.

###  3.2 Cross-correlation function

The RV is estimated as the horizontal distance between an observed mask and
the template mask of a star without a companion planet orbiting (e.g., the
horizontal distance between each square and circle in Figure 10). For the CCF
procedure, it is assumed that the template mask (of the star without a
companion planet orbiting) is predetermined, and so only the observed mask
needs to be estimated for each observed spectrum; this is done using the CCF.
First, note that the squares, circles, and triangles in Figure 10 represent
theoretical masks (i.e., compositions of Dirac $\delta$ functions),
respectively, on a continuum of wavelength values. In practice, however, an
empirical mask is constructed as a composition of indicator functions (rather
Dirac $\delta$ functions), taking value 1 on intervals each the width of a
pixel and zero else. Recall that each pixel spans some range of wavelengths,
and only the central value is recorded as the pixel wavelength for an observed
spectrum of a star at some point in time. The number of support intervals for
the mask is the same as the predetermined template mask (of the star without a
companion planet orbiting), and the locations of the support intervals are a
function of the velocity shift from the template mask. See Figure 11 for an
illustration of an example mask overlaying an observed spectrum.

Figure 11: The step function is an example of an observed spectrum of a star
for a grid of pixels, and the vertical bars denote the $n=2$ support intervals
of a mask at some velocity shift $\vartheta$. The value
$\Delta_{4,1}(\vartheta)$ is the area of the shaded region corresponding to
the fraction of pixel 4 intersected with the first support interval of the
mask, and the remaining $\Delta_{i,j}(\vartheta)$ are defined similarly.

The RV for an observed spectrum is the velocity shift that aligns the support
intervals of the mask with the absorption features (i.e., local minimums) of
the observed spectrum. The objective of the CCF procedure is to identify that
velocity shift. Assuming the spectrum is observed on a grid of $m$ pixels
indexed by $i\in\\{1,\dots,m\\}$, and assuming there are $n$ support intervals
of the mask, indexed by $j\in\\{1,\dots,n\\}$, the objective is achieved by
minimizing the $CCF$ function over a range of velocity shifts $\vartheta$,
where the $CCF$ function at $\vartheta$ is defined as the correlation between
the observed spectrum and the mask shifted by $\vartheta$ (Lafarga et al.,
2020; Baranne et al., 1996):

$CCF(\vartheta):=\sum_{i=1}^{m}\sum_{j=1}^{n}f_{i}\Delta_{i,j}(\vartheta),$

where $f_{i}$ is the observed spectrum (measured in normalized flux values),
and $\Delta_{i,j}(\vartheta)$ is the fraction of pixel $i$ intersected with
the $j$-th support interval of the mask (the $\Delta_{i,j}(\vartheta)$ are
also depicted in Figure 11). An example profile of the $CCF$ evaluated over a
range of velocity shifts is shown in Figure 12.

The CCF construction is driven by the fact that an observed spectrum will have
RV approximately given by
$v_{r}=\text{argmin}\\{\vartheta\,:\,CCF(\vartheta)\\}$. This optimization
problem is usually solved by fitting a Gaussian kernel to the CCF profile
(Baranne et al., 1996) of the form,

$\Phi(\vartheta)=b-d\cdot e^{-(\vartheta-\mu)^{2}/(2\sigma^{2})},$ (9)

where $\mu$, $\sigma$, $d$, and $b$ are scalar-valued parameters to be
learned. The parameter $\mu$ is then taken to be the RV estimate. The pink
dashed line is an example of such a fitted Gaussian kernel for the CFF profile
in Figure 12. However, Prieto (2007) argues that a quadratic model performs
better than a Gaussian or cubic model in determining the final RV estimate.
Moreover, it is possible to determine an analytical approximation for the
uncertainty of the RV estimate (Prieto, 2007; Lafarga et al., 2020).

Figure 12: An example CCF profile is displayed as the solid blue line and is
overlayed with a fitted Gaussian kernel (pink dashed line). The green dots are
the calculated bisectors (i.e., average) of the CCF profile at depths ranging
from 5% to 95% of the amplitude. The length of the horizontal black line is
the FWHM.

A strength of the CCF profile is the inherent features it contains that can be
used to determine whether the observed RV was induced by a planet or by
stellar activity (Lafarga et al., 2020). The full-width half maximum (FWHM)
and the contrast are used to measure temporal changes in the CCF profiles
derived from spectra observed at different times (Lafarga et al., 2020). The
FWHM measures the width of the average absorption feature; in Figure 12, the
FWHM is the length of the horizontal black line. The contrast is a relative
measure of the depth of the absorption feature and is defined as the amplitude
of the fitted Gaussian kernel divided by its baseline level (i.e., $d/b$ in
equation (9)).

Traditionally, it has been assumed that the RV is induced by stellar activity
if the CCF-inferred RV is correlated with other features of the CCF profile
(Simola et al., 2019; Lafarga et al., 2020). More recently, Collier Cameron et
al. (2021) have shown that Doppler shifts can be disentangled from stellar
activity using a principal component analysis of the autocorrelation function
of the CCF profile. A major weakness of this approach, though, is how to
determine the false alarm probability of a detected RV with respect to the
null hypothesis that the true RV is zero (i.e., no orbiting planet). Baranne
et al. (1996), (Prieto, 2007), and Lafarga et al. (2020) provide measurement
uncertainty quantification of the CCF parameter estimates based on Monte Carlo
experiments or a fitted model to the CCF profile. However, it is unclear how
to use this uncertainty quantification to determine statistical significance
since the distribution of the minimum value (or minimizer) of the CCF profile
is not known.

###  3.3 A linear regression method

(a)

(b)

Figure 13: (a) The difference in flux values between the blue shifted observed
spectrum and the true spectrum in Figure 10 if there was only one absorption
feature. (b) The first order Hermite Gaussian with location $\mu=4469.2$ and
scale $\sigma=0.07$ that would be fit to the difference flux in subfigure (a).

The Hermite Gaussian RV (HGRV) approach proposed by Holzer et al. (2020) is a
modern update of the CCF function. This method is motivated by the idea many
absorption features are well approximated by the kernel of a Gaussian
distribution (Gray, 2005), and the change in the normalized flux for a small
RV is approximately an additive shift in the wavelength, as described by
equation (7). The difference between a Gaussian kernel and its additive shift
results in a function that resembles Figure 13 (a), a function that is almost
identical in shape to that of the first order Hermite Gaussian function shown
in Figure 13 (b). The Hermite Gaussian function is defined as

$\psi_{1}(\lambda;\mu,\sigma):=\frac{\sqrt{2}(\lambda-\mu)}{\sigma^{3/2}\pi^{1/4}}e^{\frac{-(\lambda-\mu)^{2}}{2\sigma^{2}}},$

for some scalar-valued location and scale parameters, $\mu$ and $\sigma$,
respectively (Holzer et al., 2020). Holzer et al. (2020) formalizes this
relationship in the following steps: (1) estimate the spectrum of a star
without a companion planet orbiting (referred to as the template spectrum, but
not to be confused with the template mask); (2) identify the absorption
features in the template spectrum; (3) fit a Gaussian kernel to each
absorption feature $j\in\\{1,\dots,n\\}$; (4) let $y_{i}$ (called the
difference flux) be the difference between the flux at pixel
$i\in\\{1,\dots,m\\}$ of the observed spectrum of the star (at some point in
time) and the template spectrum; and (5) fit a simple linear regression of the
$n$ first order Hermite Gaussian curves from the fitted Gaussian kernels to
the response points $y_{1},\dots,y_{m}$, as in the model formulation:

$\begin{split}y_{i}&=v_{r}\sum_{j=1}^{n}\frac{\pi^{1/4}d_{j}\mu_{j}}{c\sqrt{2\sigma_{j}}}\psi_{1}(\lambda_{i};\mu_{j},\sigma_{j})+\epsilon_{i},\\\
\end{split}$ (10)

where $\lambda_{i}$ is the observed wavelength at pixel $i\in\\{1,\dots,m\\}$,
$\mu_{j}$, $\sigma_{j}$, and $d_{j}$ are defined as in equation (9), $c$ is
the speed of light, and $\epsilon_{1},\dots,\epsilon_{m}$ are assumed random
perturbations.

Figure 14: Illustration of Gaussian kernel (blue line) fitted to a synthetic
spectrm (black dots) and corresponding absorption features (identified by the
red horizontal lines).

In Holzer et al. (2020), the template spectrum of a star is estimated by
overlaying all of the normalized, barycentric-corrected, observed spectra and
fitting a local quadratic regression curve with weights determined from a
Gaussian kernel. The regression curve is fit to 8Å of the combined spectra at
a time, and the bandwidth parameter is determined via a generalized cross-
validation procedure. An example of the fitted template spectrum to synthetic
data is shown in Figure 14. Additionally, (Holzer et al., 2020) provides
algorithms for identifying absorption features and fitting the Gaussian kernel
parameters to each absorption feature.

There are two main advantages of the HGRV method compared to the CCF method.
First, the template spectrum used in the HGRV approach to estimate the RV is
specific to the star of interest which takes out much of the guesswork
involved in defining a mask. Second, because the RV in the HGRV approach is
estimated from a simple linear regression formulation, a false alarm
probability for the detection of a new exoplanet is well-defined and readily
calculated as the p-value on the $v_{r}$ coefficient in equation (10).

Still, there are a few major drawbacks to the HGRV method. First, the template
spectrum is assumed to estimate the spectrum of a star with no stellar
activity, which is unrealistic. Second, absorption features may not be
Gaussian shaped due to stellar activity (Simola et al., 2019; Rajpaul et al.,
2020) or significant blending; Holzer et al. (2020) attempts to account for
blending by fitting a sum of two Gaussian kernels. However, if neither model
fits well, then the absorption feature is discarded from the model. While
Holzer et al. (2020) examines this issue briefly, more exploration is needed
and improvements can be made to account for the lack of Gaussianity, perhaps
by employing a skew-normal density to fit to the absorption features.

###  3.4 Gaussian process methods

Though the HGRV method is designed to estimate RV values below 1 m/s, it
implicitly assumes stellar activity will be absorbed into the template, at
best, and, at worst, completely ignores how stellar activity can effect
downstream estimation. The application of the GP framework has gained traction
in recent years to specifically address this issue.

Rajpaul et al. (2015) and Czekala et al. (2017) independently proposed using
GPs to disentangle a corrupted RV into stellar activity and Doppler shift
components. Rajpaul et al. (2015) suggests using GPs to model the stellar
activity component of the RV, the bisector inverse slope, and another stellar
activity indicator as the dependent variables of a single time-varying
process. Jones et al. (2017) extend the work of Rajpaul et al. (2015) by
introducing new stellar activity indicators that are orthogonal to the
observed Doppler shift, extending the framework to a larger class of models,
and proposing a model selection procedure. Alternatively, Czekala et al.
(2017) introduces a template-free approach to separate the observed spectra of
a binary-star system. In fact, Rajpaul et al. (2020) extend the template-free
approach of Czekala et al. (2017) to estimate relative RVs between observed
spectra of the same star.

While these methods show promise, they all rely on a univariate regression
formulation of a multivariate regression problem, by stacking the response
variables of interest into a single response vector. Since the underlying
assumption is that all of the response variables are outputs of the same time-
varying process, it seems more appropriate to estimate these responses via
multivariate linear regression; such an approach has been developed for GPs by
Wang & Chen (2015) but has not yet been employed for exoplanet detection.

###  3.5 Overarching limitations of radial velocity methods

In addition to those already discussed, there are two overarching limitations
of radial velocity methods. First, none of the methods can be easily scaled to
handle large amounts of data since the CCF method depends on manually choosing
a mask and the HGRV and GP methods are computationally expensive. This is of
particular concern given the large amount of data expected in the coming
years. Second, the most popular methods remove tellurics from observed spectra
which removes information needed for estimating extremely small RVs (Bedell et
al., 2019; Rajpaul et al., 2020). This is of particular concern for
disentangling stellar activity noise from extremely-small Doppler shifts.

###  3.6 Empirical study of radial velocity methods

In this section, we provide a small simulation study to illustrate the
relative performances of the HGRV and CCF methods using data from Jones et al.
(2017). Jones et al. (2017) use the open-source SOAP 2.0 program (Dumusque et
al., 2014) to construct a realistic dataset of the Sun with a single sunspot
observed at different rotational phases. These simulated spectra act as a
ground-truth dataset which would be unavailable for any other star (Dumusque
et al., 2014; Jones et al., 2017), and thus, these spectra are particularly
useful for comparing the properties of the HGRV and CCF methods.

Figure 15: Monte Carlo error measurements from applying the HGRV and CCF
methods to the Doppler shifted spectra simulated from the SOAP 2.0 output. The
top, middle, and bottom panels show the root mean square error, bias, and
standard deviation, respectively. The vertical axis is the K value used for
the RVs and is the maximum of the possible observed RV (absolute value). The
SNR ratio of the spectra is on the horizontal axis. Each pair of SNR and K
values involved 100 independent simulations.

We use a subset containing 805 sampled wavelengths ranging from 4468 Å to 4476
Å of the SOAP 2.0 spectra. For each of the sampled wavelengths, there are 126
observations made at different rotational phases. Using these observations, we
build a template spectrum using the same technique as that for the HGRV
method. Next, we Doppler shift the template spectrum at a particular time $t$
by $(1-v_{r}(t)/c)\lambda)$ where $c$ is the speed of light and $\lambda$ is
the wavelength.

Since we are primarily interested in evaluating the methods when the data are
generated under small RVs and different signal-to-noise ratios (SNRs), we
arbitrarily set the orbital parameters to be equal to those of Earth and
varied $K$ to obtain different maximum values of the induced RV; specifically
in equation (8), we set $\zeta$ to be 0.0167086, $\omega$ to be 288.1, and
$\gamma$ to be 0.

Additionally, $\phi(t)$ is randomly drawn from a uniform distribution between
0 and 2$\pi$ to obtain 30 time samples, and $K$ is set to one of values in the
set $\\{0.01,0.03,0.05,0.07,0.1,1,10\\}$. The number of observed spectra was
set to 30 as Holzer et al. (2020) showed a highly accurate template could be
reconstructed from that many observed spectra. The template is then Doppler
shifted for each value of $v_{r}(t)$ calculated using these parameters to
obtain the simulated observed spectra. Following Holzer et al. (2020), we
explore SNRs ranging from 100 to 300 in steps of 50 where
$SNR:=\frac{E(f_{t})}{\sqrt{V(f_{t})}}$ (Priestley, 1981). Finally, the CCF
and HGRV methods are employed as described in Sections 3.2 and 3.3, with the
exception that the CCF mask is constructed simply as an indicator for whether
an absorption feature mean was observed at a given pixel or not. This
difference results in the absorption features not being exactly centered as
Lafarga et al. (2020) prescribed, but the difference is less than 2% on
average.

For each combination of $K$ and SNR, 100 synthetic data sets are generated,
and correlated Gaussian noise is added to the Doppler shifted spectra. We set
the $(i,j)$-th entry of the Gaussian noise covariance matrix to be

$a^{2}\exp\bigg{(}\frac{-c^{2}}{8l^{2}}\left\lvert\frac{\lambda_{i}-\lambda_{j}}{\lambda_{i}+\lambda_{j}}\right\rvert^{2}\bigg{)},$
(11)

where $c$ is the speed of light, $a$ is the amplitude set to $\bar{f}/SNR$
where $\bar{f}$ is the sample flux mean across all wavelengths and observation
times for a given $K$ value, $l$ sets the length scale, and $\lambda_{i}$ and
$\lambda_{j}$ are observed wavelength values; this covariance kernel is
adopted from Czekala et al. (2017). A value of $l=300$ was determined by trial
and error such that every neighboring five pixels exhibit sufficient
correlation to distinguish the data from an independent noise scenario.
Examining the Monte Carlo root mean square error, bias, and standard deviation
estimates for the CCF and HGRV methods in Figure 15, it is observed that the
HGRV method is far more precise than the CCF method for small RVs, across all
SNRs.

##  4 Concluding remarks

While transit and RV methods have been effective in detecting several thousand
exoplanets, there are several limitations. Transit methods suffer from
irregular observation times and observations with correlated noise. We find
evidence that the AOV, PDM, and Lomb-Scargle statistics are robust to random
missingness in the data, but only the Lomb-Scargle statistic is robust to
periodic missingness. Further, we find evidence that all of the statistics are
robust to non-sinusoidal signals, but none are robust to observations with
correlated noise. Finally, as illustrated in Section 2.4, the PDM and AOV
statistics may not follow the central beta and F distributions, respectively,
that are commonly used for calculating p-values, when the data are phase-
folded under an incorrect trial period. This is highly problematic because the
likelihood of false discovery of an exoplanet may be increased if the p-values
used for testing are not calibrated to the true type 1 error probability.

On the RV side, none of the methods can be easily scaled to handle the large
amounts of data expected to come over the next few years. The HGRV method is
promising but does not account for the effects of stellar activity which is
particularly important for detecting extremely small radial velocities.
Moreover, most methods remove tellurics from observed spectra which removes
information needed for estimating extremely small RVs (Bedell et al., 2019;
Rajpaul et al., 2020). More work is needed to critically evaluate the
statistical theory behind these models, to expand these models for multiplanet
systems, and to develop methods for detecting ever smaller RVs.

## References

  * (1)
  * Baluev (2008) Baluev, R. V. (2008), ‘Assessing the statistical significance of periodogram peaks’, Monthly Notices of the Royal Astronomical Society 385(3), 1279–1285.
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., Adrianzyk, G., Knispel, G., Kohler, D., Lacroix, D., Meunier, J.-P., Rimbaud, G. & Vin, A. (1996), ‘ELODIE: A spectrograph for accurate radial velocity measurements’, Astronomy and Astrophysics Supplement Series 119(2), 373–390.
  * Bedell et al. (2019) Bedell, M., Hogg, D. W., Foreman-Mackey, D., Montet, B. T. & Luger, R. (2019), ‘WOBBLE: a data-driven analysis technique for time-series stellar spectra’, The Astronomical Journal 158(4), 164.
  * Burke Jr et al. (1970) Burke Jr, E. W., Rolland, W. W. & Boy, W. R. (1970), ‘A photoelectric study of magnetic variable stars’, Journal of the Royal Astronomical Society of Canada 64, 353.
  * Collier Cameron et al. (2021) Collier Cameron, A., Ford, E., Shahaf, S., Aigrain, S., Dumusque, X., Haywood, R., Mortier, A., Phillips, D., Buchhave, L., Cecconi, M. et al. (2021), ‘Separating planetary reflex doppler shifts from stellar variability in the wavelength domain’, Monthly Notices of the Royal Astronomical Society 505(2), 1699–1717.
  * Czekala et al. (2017) Czekala, I., Mandel, K. S., Andrews, S. M., Dittmann, J. A., Ghosh, S. K., Montet, B. T. & Newton, E. R. (2017), ‘Disentangling time-series spectra with Gaussian processes: applications to radial velocity analysis’, The Astrophysical Journal 840(1), 49.
  * Davies (1977) Davies, R. B. (1977), ‘Hypothesis testing when a nuisance parameter is present only under the alternative’, Biometrika 64(2), 247–254.
  * Davies (1987) Davies, R. B. (1987), ‘Hypothesis testing when a nuisance parameter is present only under the alternative’, Biometrika 74(1), 33–43.
  * Davies (2002) Davies, R. B. (2002), ‘Hypothesis testing when a nuisance parameter is present only under the alternative: linear model case’, Biometrika pp. 484–489.
  * Dumusque et al. (2014) Dumusque, X., Boisse, I. & Santos, N. (2014), ‘SOAP 2.0: a tool to estimate the photometric and radial velocity variations induced by stellar spots and plages’, The Astrophysical Journal 796(2), 132.
  * Dworetsky (1983) Dworetsky, M. (1983), ‘A period-finding method for sparse randomly spaced observations or “How long is a piece of string?”’, Monthly Notices of the Royal Astronomical Society 203(4), 917–924.
  * Eyer & Bartholdi (1999) Eyer, L. & Bartholdi, P. (1999), ‘Variable stars: which Nyquist frequency?’, Astronomy and Astrophysics Supplement Series 135(1), 1–3.
  * Feigelson et al. (2021) Feigelson, E. D., De Souza, R. S., Ishida, E. E. & Babu, G. J. (2021), ‘Twenty-first-century statistical and computational challenges in astrophysics’, Annual Review of Statistics and Its Application 8, 493–517.
  * Gray (2005) Gray, D. F. (2005), The observation and analysis of stellar photospheres, Cambridge University Press.
  * Heck et al. (1985) Heck, A., Manfroid, J. & Mersch, G. (1985), ‘On period determination methods’, Astronomy and Astrophysics Supplement Series 59, 63–72.
  * Holzer et al. (2020) Holzer, P., Cisewski-Kehe, J., Fischer, D. & Zhao, L. (2020), ‘A Hermite-Gaussian based radial velocity estimation method’, arXiv preprint arXiv:2005.14083 .
  * Jones et al. (2017) Jones, D. E., Stenning, D. C., Ford, E. B., Wolpert, R. L., Loredo, T. J., Gilbertson, C. & Dumusque, X. (2017), ‘Improving exoplanet detection power: Multivariate Gaussian process models for stellar activity’, arXiv preprint arXiv:1711.01318 .
  * Khan et al. (2017) Khan, M. S., Jenkins, J. & Yoma, N. B. (2017), ‘Discovering new worlds: A review of signal processing methods for detecting exoplanets from astronomical radial velocity data [applications corner]’, IEEE Signal Processing Magazine 34(1), 104–115.
  * Lafarga et al. (2020) Lafarga, M., Ribas, I., Lovis, C., Perger, M., Zechmeister, M., Bauer, F., Kürster, M., Cortés-Contreras, M., Morales, J., Herrero, E. et al. (2020), ‘The CARMENES search for exoplanets around M dwarfs-radial velocities and activity indicators from cross-correlation functions with weighted binary masks’, Astronomy & Astrophysics 636, A36.
  * Lafler & Kinman (1965) Lafler, J. & Kinman, T. (1965), ‘An RR Lyrae star survey with the Lick 20-inch astrograph II. The calculation of RR Lyrae periods by electronic computer.’, The Astrophysical Journal Supplement Series 11, 216.
  * Leis (2011) Leis, J. W. (2011), Digital signal processing using MATLAB for students and researchers, John Wiley & Sons.
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., Gully-Santiago, M., Saunders, N., Cody, A. M., Barclay, T., Hall, O., Sagear, S., Turtelboom, E., Zhang, J., Tzanidakis, A., Mighell, K., Coughlin, J., Bell, K., Berta-Thompson, Z., Williams, P., Dotson, J. & Barentsen, G. (2018), ‘Lightkurve: Kepler and TESS time series analysis in Python’, Astrophysics Source Code Library.
  * Linnell Nemec & Nemec (1985) Linnell Nemec, A. & Nemec, J. M. (1985), ‘A test of significance for periods derived using phase-dispersion-minimization techniques’, The Astronomical Journal 90, 2317–2320.
  * Lomb (1976) Lomb, N. R. (1976), ‘Least-squares frequency analysis of unequally spaced data’, Astrophysics and space science 39(2), 447–462.
  * Mayor & Queloz (1995) Mayor, M. & Queloz, D. (1995), ‘A Jupiter-mass companion to a solar-type star’, Nature 378(6555), 355–359.
  * Priestley (1981) Priestley, M. B. (1981), Spectral analysis and time series: probability and mathematical statistics, Academic Press.
  * Prieto (2007) Prieto, C. A. (2007), ‘Velocities from cross-correlation: A guide for self-improvement’, The Astronomical Journal 134(5), 1843.
  * Rajpaul et al. (2015) Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S. & Roberts, S. (2015), ‘A Gaussian process framework for modelling stellar activity signals in radial velocity data’, Monthly Notices of the Royal Astronomical Society 452(3), 2269–2291.
  * Rajpaul et al. (2020) Rajpaul, V. M., Aigrain, S. & Buchhave, L. A. (2020), ‘A robust, template-free approach to precise radial velocity extraction’, Monthly Notices of the Royal Astronomical Society 492(3), 3960–3983.
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., Droettboom, M., Bray, E., Aldcroft, T., Davis, M., Ginsburg, A., Price-Whelan, A. M., Kerzendorf, W. E. et al. (2013), ‘AstroPy: A community Python package for astronomy’, Astronomy & Astrophysics 558, A33.
  * Scargle (1982) Scargle, J. D. (1982), ‘Studies in astronomical time series analysis. II-Statistical aspects of spectral analysis of unevenly spaced data’, The Astrophysical Journal 263, 835–853.
  * Schwarzenberg-Czerny (1989) Schwarzenberg-Czerny, A. (1989), ‘On the advantage of using analysis of variance for period search’, Monthly Notices of the Royal Astronomical Society 241(2), 153–165.
  * Schwarzenberg-Czerny (1997) Schwarzenberg-Czerny, A. (1997), ‘The correct probability distribution for the phase dispersion minimization periodogram’, The Astrophysical Journal 489(2), 941.
  * Simola et al. (2019) Simola, U., Dumusque, X. & Cisewski-Kehe, J. (2019), ‘Measuring precise radial velocities and cross-correlation function line-profile variations using a skew normal density’, Astronomy & Astrophysics 622, A131.
  * Stellingwerf (1978) Stellingwerf, R. F. (1978), ‘Period determination using phase dispersion minimization’, The Astrophysical Journal 224, 953–960.
  * VanderPlas (2018) VanderPlas, J. T. (2018), ‘Understanding the Lomb–Scargle periodogram’, The Astrophysical Journal Supplement Series 236(1), 16.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J. et al. (2020), ‘SciPy 1.0: fundamental algorithms for scientific computing in Python’, Nature methods 17(3), 261–272.
  * Wang & Chen (2015) Wang, B. & Chen, T. (2015), ‘Gaussian process regression with multiple response variables’, Chemometrics and Intelligent Laboratory Systems 142, 159–165.