#  Periodicity significance testing with null-signal templates: reassessment
of PTF’s SMBH binary candidates

Jakob Robnik,1 Adrian E. Bayer,2,3 Maria Charisi,4,5 Zoltán Haiman,6,7 Allison
Lin,6 Uroš Seljak1,8
1Department of Physics, University of California, Berkeley, CA 94720, USA
2Department of Astrophysical Sciences, Peyton Hall, Princeton University,
Princeton, NJ 08544, USA
3Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York NY 10010, USA
4Department of Physics and Astronomy, Washington State University, Pullman, WA
99163, USA
5Institute of Astrophysics, FORTH, GR-71110, Heraklion, Greece
6Department of Astronomy, Columbia University, New York, NY 10027, USA
7Department of Physics, Columbia University, New York, NY 10027, USA
8Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 93720,
USA E-mail<EMAIL_ADDRESS>
###### Abstract

Periodograms are widely employed for identifying periodicity in time series
data, yet they often struggle to accurately quantify the statistical
significance of detected periodic signals when the data complexity precludes
reliable simulations (e.g., in the presence of the correlated and/or non-
Gaussian noise). We develop a data-driven approach to address this challenge
by introducing a null-signal template (NST). The NST is created by carefully
randomizing the period of each cycle in the periodogram template, rendering it
non-periodic. It has the same frequentist properties as a periodic signal
template regardless of the noise probability distribution, and we show with
simulations that the distribution of false positives is the same as with the
original periodic template, regardless of the underlying data. Thus,
performing a periodicity search with the NST acts as an effective simulation
of the null (no-signal) hypothesis, without having to simulate the noise
properties of the data. We apply the NST method to the supermassive black hole
binaries (SMBHB) search in the Palomar Transient Factory (PTF), where Charisi
et al. had previously proposed 33 high signal to (white) noise candidates
utilizing simulations to quantify their significance. Our approach reveals
that these simulations do not capture the complexity of the real data. There
are no statistically significant periodic signal detections above the non-
periodic background. To improve the search sensitivity we introduce a Gaussian
quadrature based algorithm for the Bayes Factor with correlated noise as a
test statistic, in contrast to the standard signal to white noise. We show the
method is accurate and much faster than MCMC, which enables the computation of
the Bayes Factor for all the quasars in the PTF sample. We show with
simulations that this improves sensitivity to true signals by more than an
order of magnitude. However, using the Bayes Factor approach also results in
no statistically significant detections in the PTF data.

###### keywords:

quasars: supermassive black holes, methods: statistical, software: data
analysis

††pubyear: 2020††pagerange: Periodicity significance testing with null-signal
templates: reassessment of PTF’s SMBH binary candidates–References

##  1 Introduction

Detecting periodicity in time series data $X$ is an important task in many
scientific applications, ranging from astronomy to economics. In the
frequentist approach to hypothesis testing, one constructs a test statistic
$q(X)$, which is designed to be large only when a periodic signal is present
in the data. For example, for irregularly sampled time series, the Lomb-
Scargle periodogram score is a common test statistic (Lomb, 1976; Scargle,
1982; VanderPlas, 2018). For white noise containing a sinusoidal signal, the
Lomb-Scargle score equals the likelihood ratio between the signal and no-
signal hypotheses. However, for more complex noise and signals the
correspondence to the likelihood ratio no longer holds, rendering it a sub-
optimal test statistic.

A test statistic is calibrated by its distribution under the null hypothesis,
$\mathcal{H}_{0}$, i.e. the hypothesis that there is no signal. A signal is
detected if the p-value = $P(q>q(X)|\mathcal{H}_{0})$ is sufficiently small.
Typically, multiple simulations of the data under the null hypothesis are
generated and the p-value is estimated as the fraction of simulations for
which $q>q(X)$. This approach is robust to the assumptions about the
alternative hypothesis, i.e. hypothesis that the signal is present, and the
choice of the test statistic: sub-optimal choices result in a lower
probability of detecting the signal, but one still can be confident about the
false positive probability. However, the approach is not robust to the
assumptions about the null hypothesis. Realistic datasets are often too
complex to be accurately simulated, often leading to an underestimation of the
p-value, which in turn leads to false detections. The situation is even worse
if one attempts to reduce the computational complexity and replace simulations
with analytic estimates. For example, the analytic expressions provided in
Baluev (2008, 2013) are only valid for independent, identically distributed
(i.i.d.) Gaussian noise. Bayer & Seljak (2020) and Robnik & Seljak (2022)
proposed to estimate the p-value from the Bayes Factor, but this is also
sensitive to assumptions of the underlying model, since the Bayes Factor
depends on it.

Data-driven approaches to quantify the p-value help to avoid this issue. For
example, Bayer et al. (2021) estimates the p-value of the heighest peak in the
periodogram using the secondary peaks in the periodogram. An alternative is to
directly use the data to create effective null “simulations”. One such
approach is bootstraping (Ivezic et al., 2019a), where each permutation of the
measurements is treated as a different realization of the null. However, this
approach is only valid if the null is composed of independent and identically
distributed data.

In this work, we propose a novel data-driven approach for quantifying the
false positive probability in periodicity searches based on the modification
of the periodic search template, rather than the data. The null-signal
template (NST) is constructed from the periodic template by stretching and
compressing the individual cycles of the signal. NST will assign low
significance to the real (periodic) signal, but will have the same statistical
properties as the original template on the null. Analysis with the null-signal
template can therefore be used as an effective null simulation. This approach
has been used in other fields. For example, gravitational waves trigger a
coincident detection across multiple detectors, with a small delay, equal to
the light crossing time between the detectors. Therefore, by fixing a large
non-physical delay between the detectors, one can be sure to find only false
positives (LIGO Scientific Collaboration et al., 2004). Their distribution is
the same as the false positive distribution in the original search with
coincident detection, and thus the search with the modified delay can be used
as a substitute for the null hypothesis simulations. When applying the method
to applications such as periodograms the challenge is to find NST that
reproduces the same false positive rate as the actual search, while also being
sufficiently orthogonal to the presence of true signals such that the
resulting false positive rate is only weakly affected by the possible
existence of true signal in the data.

We apply the method to the search for quasar periodicity, which can signify
the presence of sub-parsec super-massive black hole binaries (SMBHBs). In
recent years, systematic searches in large samples of quasar light curves from
time-domain surveys have revealed around 150 candidates. For instance, Graham
et al. (2015) identified 111 candidates in a sample of around 250,000 quasars
from the Catalina Real-time Transient Survey (CRTS), while Charisi et al.
(2016) identified 33 candidates among around 35,000 quasars from the Palomar
Transient Factory (PTF). More recent searches in the Panoramic Survey
Telescope and Rapid Response System (Pan-STARRS), the Dark Energy Survey (DES)
and the Zwicky Transient Facility (ZTF) detected a few more candidate systems
(Liu et al., 2019; Chen et al., 2020, 2024). For a detailed summary of the
current status of quasar periodicity searches, we refer the reader to Charisi
et al. (2022); D’Orazio & Charisi (2023) Here we apply our new method to
reassess the periodicity significance of the 33 SMBHB candidates identified in
PTF.

In addition to the NST method, which is agnostic to the underlying noise
properties, we also introduce some improvements to the standard quasar
periodicity detection methods, which typically model quasar variability as a
damped random walk (DRW) (MacLeod et al., 2010; Kelly et al., 2009; Kozłowski
et al., 2010). Previous searches for quasar periodicity (Graham et al., 2015;
Charisi et al., 2016; Chen et al., 2020) used periodogram peaks as their test
statistic and then assessed the p-values of the detected signals with DRW
simulations. Given that the underlying assumption in the periodogram is
uncorrelated white noise, this makes the choice of this statistic sub-optimal,
as the noise is correlated. In addition, priors may also improve the search
sensitivity, by focusing on the regions of parameter space with higher prior.
The posterior odds, often equivalent to the Bayes Factor, are known to be the
optimal test statistic (Zhang, 2017; Fowlie, 2023) in the sense that they
maximize the true positive rate (TPR) at a fixed false positive rate (FPR).
However, the Bayes Factor with correlated Gaussian noise has been computed
only for a few particularly interesting candidates (D’Orazio et al., 2015; Zhu
& Thrane, 2020; Vaughan et al., 2016). This is because the evidence integral
is usually computed with sampling algorithms such as nested sampling (Zhu &
Thrane, 2020; Witt et al., 2022), which is computationally intensive.
Sometimes a Bayesian information criterion (BIC) is used as an easy to
evaluate substitute for the Bayes Factor (Witt et al., 2022; Liu et al.,
2018), but BIC does not take into account the prior information and accounts
for the look elsewhere effect overly simplistically leading to a suboptimal
analysis (Weakliem, 1999).

Exploiting the low dimensionality of the problem, we instead compute the Bayes
Factor with a Gaussian quadrature scheme (Robnik & Seljak, 2022), achieving
the same accuracy as sampling methods at about three orders of magnitude lower
computational cost. Specifically, a single light curve analysis with our code
takes around half a minute on a single CPU, while in Witt et al. (2022),
several CPU hours were used. The speed-up enables the use of the Bayes Factor
to search for periodicity in the entire quasar sample, which previously would
have been prohibitive due to computational demands. This can be scaled to the
even larger samples in future surveys, like the Rubin Observatory (Ivezic et
al., 2019b; Xin & Haiman, 2021) or the Roman telescope (Haiman et al., 2023),
which will detect millions of quaasars.

In Section 2 we present the methodology of the novel data-driven approach to
use null-signal templates to create effective null simulations. In Section 3
we validate the method on synthetic light curves. In Section 4 we then present
the improvements to the SMBHB search in PTF data and apply the NST to
calibrate the results. The code with tutorials is available at
https://github.com/JakobRobnik/periodax.

##  2 Method

Let us consider a time-series $X$, consisting of measurements
$\\{x_{i}\\}_{i=1}^{N}$, with measurement errors $\\{\sigma_{i}\\}_{i=1}^{N}$,
taken at times $\\{t_{i}\\}_{i=1}^{N}$, where $N$ is the total number of
observations. Let $T$ be the time span of the data (i.e the baseline of the
light curve). The measurements $x_{i}$ are composed of noise $n_{i}$ and
signal $s(t_{i})$, so that:

$x_{i}=s(t_{i})+n_{i},$ (1)

where the signal

$s(t)=A\,u\left(\frac{t}{P}-\phi\right)$ (2)

is periodic, i.e. $u(x+k)=u(x)$ for integer $k$. $P$ is the period,
$\phi\in[0,1)$ the phase, and $A$ the amplitude of the periodic signal. The
standard periodogram corresponds to a simple sinusoidal signal,
$u(x)=\sin\\{2\pi x\\}$. Other examples include exoplanet transit searches,
where $u(x)$ is a U-shaped transit (Robnik & Seljak, 2021), or more complex
SMBHB signals, like sawtooth periodicity (Duffell et al., 2020; Westernacher-
Schneider et al., 2022) or periodically repeating lensing flares (D’Orazio &
Stefano, 2020).

We would like to know if a periodic signal is present in the data. The
strategy is to design some test statistic $q(X)$. A popular test statistic for
periodicity detection is the Lomb-Scargle periodogram score:

$\displaystyle q_{LS}(X|s)$
$\displaystyle=\max_{A,P,\phi}\,\sum_{i,j=1}^{N}x_{i}\Sigma^{-1}_{ij}x_{j}-\Delta_{i}\Sigma^{-1}_{ij}\Delta_{j}$
(3) $\displaystyle\Delta_{i}\equiv x_{i}-s(t_{i}|A,P,\phi),$

where $\Sigma$ is the noise covariance matrix, which is diagonal in the
standard white Lomb-Scargle periodogram,
$\Sigma_{ij}=\delta_{ij}\sigma_{i}^{2}$, and only accounts for the measurement
errors. Another common test statistic is the likelihood ratio

$q_{LR}=2\log\frac{\max_{\boldsymbol{z}_{1}}p(X|\mathcal{H}_{1},\boldsymbol{z}_{1})}{\max_{\boldsymbol{z}_{0}}p(X|\mathcal{H}_{0},\boldsymbol{z}_{0})},$
(4)

where $\boldsymbol{z}_{i}$ are the parameters of the $\mathcal{H}_{i}$
hypothesis. Lomb-Scargle score is a special case of the likelihood ratio,
namely when the noise is Gaussian and the null has no parameters. While the
likelihood ratio takes the ratio of the likelihoods at the optimal parameters,
the Bayes Factor takes the ratio of the evidences:

$B=\frac{p(X|\mathcal{H}_{1})}{p(X|\mathcal{H}_{0})},$ (5)

which are the likelihoods averaged over the prior $p(X|\mathcal{H}_{i})=\int
p(X|\boldsymbol{z}_{i})p(\boldsymbol{z})d\boldsymbol{z}_{i}$. Finally, the
posterior odds additionally weigh the Bayes Factor with the prior odds, i.e.
the ratio $p(\mathcal{H}_{1})/p(\mathcal{H}_{0})$ that favors the datasets
where the signal is more probable a priori:

$\frac{p(\mathcal{H}_{1}|X)}{p(\mathcal{H}_{0}|X)}=\frac{p(X|\mathcal{H}_{1})}{p(X|\mathcal{H}_{0})}\frac{p(\mathcal{H}_{1})}{p(\mathcal{H}_{0})}.$
(6)

Once a test statistic is chosen, it is evaluated on the data $q(X)$ and is
compared to the distribution $P(q|\mathcal{H}_{0})$ under the null hypothesis.
If the observed value is very unlikely under the null, we can claim a
discovery.

Here we propose to estimate the null distribution by modifying the test
statistic $q(X)\xrightarrow[]{}\widetilde{q}(X)$, such that the statistical
properties of both statistics are the same under the null:

$P(\widetilde{q}|\mathcal{H}_{0})\approx P(q|\mathcal{H}_{0}),$ (C1)

while the modifed statistic is only weakly triggered on the alternative

$P(\widetilde{q}|\mathcal{H}_{1})\approx P(\widetilde{q}|\mathcal{H}_{0})\ll
P(q|\mathcal{H}_{1}).$ (C2)

A test statistic satisfying Conditions (C1) and (C2) can be used as an
effective null simulation when applied to the real data, regardless of whether
the data is actually from $\mathcal{H}_{0}$ or $\mathcal{H}_{1}$. This is a
purely data driven approach, and does not require any simulations to quantify
the false positive rate.

Figure 1: An example of a sinusoidal signal template (teal) and the null
sinusoidal template (orange) from Eq. (7). The individual cycles of the null-
signal template are marked with arrows.

We construct $\widetilde{q}$ by perturbing the signal template of the
alternative hypothesis $\mathcal{H}_{1}$, such that
$\widetilde{q}(X)=q(X|\widetilde{s})$. For example, $\widetilde{q}_{LR}(X)$ is
the ratio of the likelihoods, but the alternative hypothesis uses the
perturbed template. The null-signal template $\widetilde{s}$ is the same as
the original template on short time scales, i.e. the period of the signal, but
differs significantly on time scales comparable to the length of the light
curve. The similarity of the templates on short time scales ensure that $q$
and $\widetilde{q}$ trigger similarly on the null signal, satisfying condition
(C1), while long-scale differences ensure that the overlap between the null-
signal template and the original template is small, hence satisfying condition
(C2). This method thus works best when the number of signal cycles is large,
but we provide an additional prescription to apply it with as few as two
signal cycles.

There are multiple ways of modifying the template; in this work, we randomize
the length of its cycles. The original template is periodic, so the $n$-th
cycle is the interval $[(n-1)P,\,nP]$. The null-signal template’s cycle $n$
will have length $P_{n}$, such that the $n$-th cycle now spans the interval
$[\sum_{k=1}^{n-1}P_{k},\,\sum_{k=1}^{n}P_{k}]$. Inside each cycle, the null-
signal template equals the stretched or contracted original template cycle.
More formally, the null-signal template is

$\widetilde{s}(t)=A\,u\left(\frac{t-t_{0}(t)}{P_{i(t)}}-\phi\right),$ (7)

where $i(t)\equiv\mathrm{max}\\{j|\sum_{k=1}^{j}P_{k}<t\\}$ and
$t_{0}(t)\equiv\sum_{k=1}^{i(t)}P_{k}$. Note that this template is continuous
if $u$ is continuous, but it is not differentiable. An example of the null-
signal template is shown in Figure 1 in orange, while the original periodic
template is shown in teal.

It remains to select the lengths $P_{n}$ of cycles. The original template Eq.
(2) is a function of the signal period $P$, so an equivalent of the period
must also be defined for the null-signal template Eq. (7). Furthermore, the
null-signal template must vary smoothly with the period. We achieve this by
requiring that the null-signal template at period $P$ has the same number of
possibly noninteger cycles as the original template with period $P$, which is
$T/P$. This amounts to the following constraint

$P_{1}+P_{2}+\ldots P_{n}+\lambda P_{n+1}=T=(n+\lambda)P.$ (8)

Here, the number of cycles is decomposed into integer and fractional parts:
$T/P=\lfloor T/P\rfloor+\text{frac}(T/P)\equiv n+\lambda$. The constraint from
Eq. (8) is satisfied if we pick

$P_{n}=\frac{Tp_{n}}{\lambda p_{n+1}+\sum_{k=1}^{n}p_{k}},$ (9)

where $\\{p_{n}\\}$ can be any positive numbers. For example, one could use
i.i.d. draws from the Gamma distribution which would result in periods being
symmetrically Dirichlet distributed. We would like to avoid the tails of the
Gamma distribution so we pick the uniform distribution,
$p_{n}\sim\mathcal{U}(1/\sqrt{3},\sqrt{3})$.

Figure 2:  Lomb-Scargle periodogram score with the null-signal template (NST)
relative to the score with the periodic template, i.e. fractional overlap
$r(s_{\mathrm{NST}},s_{\mathrm{periodic}})$. The data here are just the
periodic template with no noise. Lower is better, because it ensures NST has
small overlap with the periodic template and consequently obeys Condition
(C2). The red line corresponds to NST if all periods are drawn from the
uniform distribution. The blue line uses the fixed first seven periods that
minimize the overlap and random other periods, as described in Section 2. Time
sampling is 20 equally spaced points / cycle and the frequency grid for the
NST is 1000 logarithmically spaced frequencies in the range $[P/10,10P]$. The
solid lines are medians over two hundered NST period randomness realizations.
The stronger and the weaker band correspond to 25-75 and 10-90 quantiles
respectively. As can be seen the optimized periods (blue) correspond to the
best possible case of the random periods (red) and reduce the overlap to
$60\%$ even at only two signal cycles.

Having random periods can be problematic if the number of cycles is very low,
because all periods can become very similar to each other by pure chance.
Should this happen, the null-signal template would no longer have small
overlap with the periodic template and would violate Condition (C2). To
prevent this from happening we choose $p_{n}$ deterministically for low $n$.
We define the fractional overlap between two templates as
$r(s_{\mathrm{test}},s_{\mathrm{true}})=q_{LS}(s_{\mathrm{true}}|s_{\mathrm{test}})/q_{LS}(s_{\mathrm{true}}|s_{\mathrm{true}})$.
This is the ratio of the Lomb-Scargle scores that the template
$s_{\mathrm{test}}$ would assign to a signal which was actually generated by
$s_{\mathrm{true}}$, relative to the score that we would get using the
periodic template $s_{\mathrm{true}}$. Small values of the fractional overlap
will ensure that Condition (C2) holds. We can achieve this by optimizing the
values $p_{n}$. We do this iteratively, starting with $p_{1}=1$ without loss
of generality. In the $n$-th iteration, $p_{1},p_{2},\ldots p_{n}$ have
already been fixed. We take the time span $T=(n+1)P$ and optimize $p_{n+1}$
such that the overlap $r(s_{\mathrm{true}},s_{\mathrm{test}})$ is minimized.
To make the optimized values more robust to the possible gaps in the data, we
try $n$ different gap options, each covering one of signal cycle. The function
that we then minimize is $r$, maximized over the gap choice. To compute the
overlap we use a sinusoidal template with an almost continuous time sampling
(20 equally spaced points / cycle) and a dense grid for the periodic template
trial frequencies (1000 logarithmically spaced frequencies in the range
$[P/10,10P]$). We optimize for the first 7 $p_{n}$ and normalize them to unit
average to match the scale of $p_{n>7}$. These are then fixed throughout the
paper. Figure 2 shows that this procedure yields small overlaps, even when the
time span of the data is not much larger than the period of the signal.

Note that this method only ensures that the fractional overlap between the
periodic template and the NST is small, but not its absolute value. Therefore,
the presence of a very large alternative hypothesis signal can compromise the
first equality in condition (C2). In this case, it may be better to perform
iterative analysis, such as subtracting the best fit signal from the data and
compute the modified test statistic on residuals, otherwise the p-value
estimate would be conservative. Either way this case is easy to identify, and
will not be an issue for the current work, since as we will show no periodic
signal candidates surpass the NST baseline.

##  3 Validation

Figure 3: Time sampling in synthetic datasets from Section 3. Time is on the
x-axis (but with different scales for different time samplings), the dots
indicate times where measurements are taken. From the top to the bottom:
equally spaced, equally spaced with gaps, randomly spaced. The bottom row is
an example PTF time sampling for a quasar SDSS J085037.61+201337.1. Figure 4:
On noise-only simulations, the p-value with the null-signal template (colored)
agrees perfectly with the periodic template (black). A grid with all
possibilities of time sampling and noise distribution from Section 3 is shown.
This demonstrates that validity of Condition (C1) is not sensitive to the
noise and time sampling properties. Note the large case to case differences in
the signal scores. The shaded bands are statistical one $\sigma$ uncertainties
due to the finite number of simulations, obtained with bootstraping. Figure 5:
On noise with injected signal simulations, the score distribution with the
null-signal template (blue) differs significantly from the periodic signal
(orange) and matches almost perfectly the distribution with the periodic
template on noise-only simulations (black). This demonstrates that validity of
Condition (C2) is not sensitive to the noise and time sampling properties. A
grid with all possibilities of time sampling and noise distribution from
Section 3 is shown.

We now test Conditions (C1) and (C2) in a variety of synthetic examples where
exact null simulations can be performed and compared to the NST effective null
simulations. The examples are chosen to provide some variety of noise
properties and time samplings. We consider three different noise
distributions:

  * •

i.i.d. standard Gaussian,

  * •

i.i.d. Cauchy,

  * •

correlated Gaussian with a DRW kernel, $\tau=T/20$,
$\sigma=1$ and measurement errors $\sigma_{i}=0.2$. $T$ is the time span of
the light curve.

with three different time samplings:

  * •

1000 equally spaced observations, $t_{k}=k$,

  * •

equally spaced with some observations missing. We select four large regions
where the data is missing and additionally independently remove each
measurement with 10 % probability. This time sampling somewhat resembles
Kepler space telescope sampling (Jon M. Jenkins et al., 2017)

  * •

1000 uniformly random sampled $t_{k}\sim\mathcal{U}(0,1)$.

Time samplings are shown in Figure 3. For each of the 9 combinations of time
samplings and noise, we do $8192$ simulations. For the test statistic we use
the floating mean periodogram score from Appendix A, maximized over the trial
frequencies $\\{k/T\\}_{k=1}^{N/2}$. The resulting cumulative distribution
with the sinusoidal and the null template are shown in Figures 4 and 5. In
Figure 4 we perform noise-only simulations and show a perfect match between
the two templates. Thus Condition (C1) is satisfied and is not sensitive to
the noise or time sampling properties. In Figure 5 we inject a sinusoidal
signal in the data. In each realization the signal phase
$\phi\sim\mathcal{U}(0,1)$ and frequency $1/P\sim\mathcal{U}(1/T,N/2T)$ are
drawn from a uniform distribution. The amplitude of the injected signal is
chosen for each combination of time sampling and noise separately, such that
it significantly perturbs the distribution of the test statistic in Figure 4.
On the injected sinusoidal data, the null-signal template (teal) gives a
significantly lower score to the periodic template (orange) and almost
perfectly matches the pure-noise simulations with the periodic template
(black). This means that Condition (C2) is also satisfied and as before is not
sensitive to the noise or time sampling properties. We will further test these
conditions in the next section with simulations tailored to the SMBHB search
in the PTF data.

##  4 SMBHB search in PTF data

###  4.1 PTF Light Curves

We extract quasar light curves obtained by PTF (Law et al., 2009) and apply
the same cuts as in Charisi et al. (2016) to obtain a sample of 35,383
spectroscopically confirmed quasars.

Following Charisi et al. (2016), we model the null hypothesis as a damped
random walk (DRW) with additional Gaussian measurement noise and a constant
offset to represent the mean magnitude. A DRW is a correlated Gaussian noise
with covariance matrix

$\Sigma_{ij}^{(DRW)}=\sigma^{2}e^{-|t_{i}-t_{j}|/\tau},$ (10)

where $\sigma$ is the strength of the DRW noise and $\tau$ its correlation
time. We adopt a log-normal prior for both parameters,

$\log\sigma\sim\mathcal{N}(\log 0.1,0.2)\qquad\log\tau\sim\mathcal{N}(\log
120\,\text{days},0.9),$ (11)

following Figure 9 from Charisi et al. (2016) and MacLeod et al. (2010). The
total noise is a combination of the DRW process and the measurement noise,
thus the total covariance matrix is

$\Sigma_{ij}=\Sigma_{ij}^{(DRW)}+\sigma_{i}^{2}\delta_{ij}.$ (12)

We first remove the outliers as follows: we identify the data point $i$ as an
outlier if it deviates from the Gaussian Process (GP) fit to the data by more
than $3$ standard deviations,

$|x_{i}-\mu_{GP}(t_{i})|>3\sqrt{\sigma_{i}^{2}+\sigma^{2}_{GP}(t_{i})}.$ (13)

We take the GP with the kernel from Eq. (12) with fixed $\sigma=0.1$ and
$\tau=100$ days, where $\mu_{GP}$ and $\sigma^{2}_{GP}$ are the GP posterior
mean and variance, respectively. We iteratively fit the GP to the light curves
and identify the outliers for each new fit, repeating the procedure 5 times.
We then bin together the observations taken on the same night, as in Charisi
et al. (2016).

###  4.2 Signal prior

The alternative hypothesis ($\mathcal{H}_{1}$) is a sinusoidal signal with
period $P$. In the quasar reference frame the period is expected to be
distributed as $p(P_{RF})\propto P_{RF}^{\alpha}$, with $\alpha=8/3$, assuming
that SMBHBs are inspiraling on circular orbits and their evolution is driven
by gravitational waves only (Haiman et al., 2009). The observed period is
related to the reference frame period by $P=P_{RF}(1+z)$, where $z$ is the
redshift of the quasar. Therefore the observed period is also distributed as
$p(P)\propto P^{\alpha}$. We will enforce an upper limit on the period of
$P<P_{\mathrm{cut}}=T/n_{\mathrm{min}}$ to ensure that at least a minimal
number of signal cycles $n_{\mathrm{min}}=2$ are observed, otherwise
periodicity detection cannot be claimed.111Note that previous searches imposed
a more relaxed requirement of only 1.5 cycles, which may also contribute to
the detection of false positives. The probability that a given quasar has a
signal with a period in the desired range of $0<P<P_{\text{cut}}$ is

$P(\mathcal{H}_{1})\propto\int_{0}^{P_{\mathrm{cut}}/(1+z)}p(P_{RF})dP_{RF}\propto\frac{P_{\mathrm{cut}}^{\alpha+1}}{(1+z)^{\alpha+1}},$
(14)

where $P_{\mathrm{cut}}/(1+z)$ is the period cutoff translated to the quasar
reference frame. The prior odds $P(\mathcal{H}_{1})/P(\mathcal{H}_{0})$
penalizes quasars with larger redshift and shorter data span, because the
required number of observed cycles limits the allowed observed period which
translates to a smaller range of reference-frame periods. The proportionality
constant in Eq. (14) further depends on quasar’s brightness, because brighter
quasars with more massive putative binaries chirp faster, and so spend less
time at some fixed period (Peters, 1964). We will ignore this effect here and
assume the proportionality constant is the same for all quasars. We will fix
this constant to a value that ensures
$\sum_{i}P(\mathcal{H}_{1}|\mathrm{quasar\,i})=\sum_{i}P(\mathcal{H}_{0}|\mathrm{quasar\,i})$.
This is in line with the standard practice in hypothesis testing, i.e. that no
hypothesis is preferred a priori. The choice of the constant has no effect on
the posterior odds as a test statistic, because it just shifts all scores by a
constant. Finally, to make the prior smooth at the edges, we multiply $p(\log
P)$ by the cosine-tapered window which has a value of 1 for
$P<P_{\mathrm{min}}/1.2$ and 0 for $P>P_{\mathrm{min}}$. This makes it easier
to maximize the posterior density during the computation of the Bayes Factor
and also better aligns with our prior beliefs, which dictate a smooth prior.
For the phase we adopt the uniform prior $\phi\sim\mathcal{U}(0,1)$. For the
amplitude we also adopt the uniform prior, but do not specify its upper bound
and ignore the trials factor associated with the amplitude. This only amounts
to a constant shift of the test statistic which has no relevance for the Bayes
Factor as a test statistic (Robnik & Seljak, 2022).

###  4.3 Test statistics

Figure 6: True positive probability as a function of false positive
probability (ROC curve) for the periodicity search in the PTF data is shown
for the test statistics from Section 4.3. The shaded regions are $1\sigma$
confidence bands obtained by bootstrapping. The posterior odds and the Bayes
Factor perform best and are closely followed by the log-likelihood ratio.
White periodogram score performs much worse.

We try several test statistics:

  * •

Standard (uncorrelated) peridogram score, Eq. (3).

  * •

Likelihood ratio, Eq. (4)

  * •

Bayes Factor, Eq. (5)

  * •

Posterior odds, Eq. (6)

We show how these can be efficiently computed in Appendices A and B. In Figure
6 we show the receiver operating characteristic (ROC) curve for these test
statistics. For this, we generated mock light curves by taking the PTF quasar
data, described in Section 4.1 and replacing the magnitude measurements with
realizations of DRW, whose parameters are drawn from the prior, as in Eq.
(11). A signal with amplitude $A=2\sigma$ and phase and period drawn from the
priors in Section 4.2 is randomly injected in some light curves. The signal is
injected according to the prior odds, so in around half of the light curves,
but with higher probability where the prior odds are larger. To compute the
ROC curve for a given test statistic, the detection threshold is varied and at
each value the fraction of recovered true signals (i.e. true positive
probability) is plotted against the fraction of identified false signals (i.e.
false positive probability).

Figure 6 shows that the posterior odds (yellow) are optimal in the sense that
they maximize the true positive probability at a fixed false positive
probability. This is the content of the Neyman-Pearson lemma for the composite
hypothesis (Zhang, 2017; Fowlie, 2023), so we have here numerically confirmed
it. Posterior odds are almost indistinguishable from the Bayes Factor (teal),
meaning that the prior odds do not play a significant role in this
application. The likelihood ratio (orange), which accounts for the correlated
noise, but does not marginalize over the unknown parameters and does not take
into account the prior distribution, closely follows the posterior odds and
the Bayes Factor but is slightly less optimal. The standard periodogram score
(red), which assumes white noise, is significantly suboptimal: at 0.001 false
positive probability it has a 20 times lower chance of detecting a real
signal, because it needs a higher detection threshold to keep the number false
positives low. Most of the false positives arise from the correlated signal of
the quasar lightcurves, which this test statistic ignores.

The take-away message is that the detection efficiency benefits significantly
from incorporating the noise correlations in the test statistic. Because the
Bayes Factor and posterior odds perform almost identically and the Bayes
Factor is the standard choice, we will use the Bayes Factor in the remainder
of the paper. Note that one would not benefit significantly in terms of
computational cost by using the likelihood-ratio instead. This is because the
Gaussian quadrature integration is cheaper than the optimization, which is
needed to identify the optimal parameters and is needed for both test
statistics.

###  4.4 Null-signal template validation

Figure 7:  Cumulative number of detections as a function of the Bayes Factor
test statistic on a simulation of the SMBHB search in PTF data. On noise-only
simulation (DRW), the null-signal template (brown) and the periodic template
(black) match perfectly, thus Condition (C1) is satisfied. On simulations of
the noise with injected signal, the null template (orange) yields
significantly lower scores than the periodic template (red) and matches the
noise-only simulations. Thus Condition (C2) is also satisfied.

Next we repeat the test for Conditions (C1) and (C2), as in Section 3, but
specifically for the SMBHB search in PTF. For this, we take the PTF light
curves from Section 4.1 and replace the magnitude measurements with
realizations of DRW, whose parameters are drawn from the prior in Eq. (11).
These light curves serve as simulations of the null hypothesis
$\mathcal{H}_{0}$. As explained in 4.3, we then use the Bayes Factor as our
test statistic. In Figure 7 we show that the p-value computations with the
periodic signal template (black) and the null-signal template (brown) agree
perfectly, thus demonstrating Condition (C1). We also inject a sinusoidal
signal in all light curves, with properties as in the previous section, i.e.
with an amplitude of $A=2\sigma$, and phase and period drawn from their prior
(Section 4.2). This light curves serve as simulations of the alternative
hypothesis $\mathcal{H}_{1}$. Figure 7 shows that the NST distribution
(orange) produces significantly lower scores than the periodic signal (red),
and thus Condition (C2) is satisfied.

###  4.5 Results

Figure 8: The results of SMBHB search in PTF data. Cumulative count of
detections as a function of the detection threshold is shown. The periodic
template results (teal) are to be compared with the null-signal template
results (orange), which act as a null simulation. Detection does not exceed
the null baseline, so a discovery cannot be claimed. In contrast, DRW
simulations with parameters drawn from the prior (black) yield significantly
lower scores, which would lead to a false signal claim. We have performed 100
simulations of the entire dataset analysis and sorted the candidates in each
analysis by their score. The percentiles are then computed for each rank
separately.  Figure 9:  Left: The results of SMBHB search in PTF data with the
standard white Lomb-Scargle periodogram test statistic $q_{LS}$ from Equation
(3). The cumulative count of detections as a function of the detection
threshold is shown. The periodic template results (teal) are to be compared
with the null-signal template results (orange), which act as a null
simulation. Periodic template results do not exceed the null baseline in a
statistically significant manner, so a discovery cannot be claimed. The
confidence regions are the quartiles obtained by bootstrapping. Here, even the
(100) DRW simulations come very close to the null signal template results,
showing that the majority of false positives are a result of the test
statistic not incorporating the DRW noise. Right: the correlation between the
white periodogram score and the Bayes factor for all the quasars, using the
periodic template. Two test statistics are not highly correlated, especially
in the tails: high value of the Bayes factor does not imply high value of the
white periodogram score. In particular, the candidate with the highest white
periodogram score (marked in red) causes a slight deviation between the NST
and periodic template results in the left panel, but has a low value of the
Bayes factor. Thus, it is most likely caused by DRW not being properly modeled
in the white periodogram.  Figure 10: A demonstration that the null-signal
template gives a lower significance to the periodic signal. The top 30 quasars
from Figure 8 are taken. The reason why these quasars rank so highly might be
some corruption of their data, so we extract the optimal signal from these
quasars and inject them in 30 among the lower ranking quasars, which will here
act as realistic null simulations. To ensure that the signal was injected in a
similar environment, we select this second set of 30 quasars to match the
original quasars in terms of time span of the data, number of observational
nights and parameters of the DRW kernel. In 90% of cases we were able to match
all of those parameters to less than 10% deviations. We thus obtained 30
simulations of the alternative hypothesis, which is similar to the claimed
SMBHB candidates. We analyze these data sets both with the periodic template
(teal) and the null-signal template (orange). We show that NST gives lower
scores than the periodic template, so such a SMBHB population would have been
detectable.

Finally, we apply all the proposed improvements for the periodicity search, in
the PTF light curves. In Figures 8 and 9, we show the results of the
periodicity search, respectively using the Bayes Factor and white periodogram
score as the test statistic. All the 33 SMBHB candidates proposed in Charisi
et al. (2016) achieve $\log B<11$ so they are no longer among the top
candidates when the Bayes Factor is used as a test statistic. The top ranking
candidates achieve $\log B>20$. This shows that the choice of the test
statistic has a significant influence on the results. We have argued that
Bayes Factor achieves a significantly higher detection sensitivity.

Based solely on the raw Bayes Factor values, one might be inclined to confirm
the top candidates, as even after accounting for the multiplicity of trials,
$B\approx\exp{20}/35383\approx 10^{4}$, where we have treated all quasars
equally. This is much larger than $100$, a decisive evidence limit by the
Jeffrey’s interpretation scale (Jeffreys, 1998). Similarly, the p-values
obtained by DRW simulations would suggest a low p-value for these candidates,
see Figure 8. However, both of these approaches suffer from poor modeling of
the null hypothesis. This is demonstrated with the NST method, which is more
robust and agnostic to the assumptions for the underlying noise. Bayes Factors
and white periodogram scores, estimated with the NST, reveal that there are no
statistically significant periodicity detections, compared to the non-periodic
background. NST reveals a much higher Bayes Factors than DRW simulations,
demonstrating that the latter are a poor representation of the real
variability of quasars.

In fact, NST yields higher Bayes Factor scores, demonstrating that there might
be even more non-periodic than periodic signal data. In other words, non-
periodic signal would have been a better alternative hypothesis than the
periodic signal. In contrast, we show in Figure 10 that a real signal would
trigger a significantly higher detection with the periodic template than with
NST. Thus whatever the source of NST giving higher Bayes Factor than periodic
signal is, it cannot be caused by the presence of real periodic signal in the
data. It is however possible that the NST signal being above the periodic
signal is caused by the presence of quasi periodic signals in the data that
have templates similar to NST. However, investigating this hypothesis would
require a more detailed analysis which is beyond the scope of this paper,
where we only focus on periodic signals.

In Figure 11 we show 5 light curves with the highest Bayes Factor, both with
the periodic and the null-signal template. In most cases, the lightcurves are
similar, suggesting the periodogram is triggered on false positives.

##  5 Summary

We proposed a novel method to improve the accuracy of detecting periodic
signals in complex time series data by modifying the periodogram template to
create a non-periodic null-signal template (NST). Since anything detected by
the NST is by default a false positive, this allows us to use the real data as
effective null simulations, avoiding the need for potentially inaccurate data
simulations. This is particularly useful in fields like astronomy, where the
data complexity and the noise patterns are challenging to model accurately.

We validated our method through synthetic examples and applied it to the
search for quasar periodicity, which have been suggested to track the presence
of SMBHBs, in $\sim$35,000 quasar light curves from PTF. Our results showed
that NST provides a robust estimate of the false positive rate, and revealed
that previously proposed SMBHB candidates are likely false positives. The
false positive rate in previous analyses was likely underestimated due to the
assumption of DRW variability; simulations of DRW light curves cannot
reproduce the NST distribution, which indicates that the DRW model may be a
limited description of quasar variability, for example that its correlation
time may be longer than expected, or that the distribution is non-Gaussian.
Since all the systematic searches for quasar periodicity to date assume DRW
model for quasar variability (Graham et al., 2015; Liu et al., 2019; Chen et
al., 2020, 2024), many of the identified candidates in these searches may also
be false positives, and would need to be validated with the NST method. We
plan to apply our method to light curves from other surveys, like CRTS and
ZTF.

Additionally, we compared several test statistics (periodogram score, Bayes
Factor, posterior odds, and log likelihood ratio). We demonstrated that the
Bayes Factor (and posterior odds) outperform the periodogram score, enhancing
the sensitivity of periodicity searches by a factor of 20. This indicates that
including the correlated noise in the test statistic is beneficial. In the
future work the test statistic could have been improved further, by computing
the Bayes Factor against the true null hypothesis, which is not captured by
DRW, as suggested by our results.

We also introduced a faster algorithm for computing the Bayes Factor of
sinusoidal SMBHB signals against the DRW quasar noise, enabling its
application to large datasets, like the upcoming LSST of the Rubin
Observatory.

Looking towards future applications of NST, it is worth repeating the analysis
with non-sinusodial periodic templates, as well as with quasi-periodic
templates. The challenge of these applications is again to find NSTs that
satisfy conditions C1 and C2. Furthermore, while we have applied NST to a
problem concerning SMBHB, the method will be useful in any application where
one is searching for periodic signals, for example in exoplanet detection.

Figure 11: Top: 5 quasars with highest Bayes Factors for the SMBHB hypothesis.
Time series data and the best fit models are shown. Time is measured relative
to MJD 54903 days. Bottom: same, but with the null-signal template.

## Acknowledgements

This material is based upon work supported in part by the Heising-Simons
Foundation grant 2021-3282 and by the U.S. Department of Energy, Office of
Science, Office of Advanced Scientific Computing Research under Contract No.
DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory to enable research
for Data-intensive Machine Learning and Analysis. MC is funded by the European
Union (ERC-STG-2023, MMMonsters, 101117624). ZH acknowledges support by NSF
grant AST-2006176 and NASA grants 80NSSC22K0822 and 80NSSC24K0440.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

##  Appendix A Periodogram as a matched filter

The likelihood of the Gaussian data under the null hypothesis is

$-2\log
p(\boldsymbol{x}|\Sigma)=\boldsymbol{x}^{T}\Sigma^{-1}\boldsymbol{x}+\log\det
2\pi\Sigma$ (15)

and

$-2\log p(\boldsymbol{x}|\boldsymbol{z},\Sigma)=-2\log
p(\boldsymbol{x}-\boldsymbol{s}(\boldsymbol{z})|\Sigma)$ (16)

under the alternative hypothesis. $\Sigma$ is the noise covariance matrix
which we will assume to be known in this section. $\boldsymbol{z}$ are the
parameters of the signal. Let the signal depend on the first $M$ parameters
linearly, we will call these parameters the amplitudes and denote them by
$\boldsymbol{a}$. The signal template is then of the form

$\boldsymbol{s}=\sum_{i=1}^{M}a^{i}\boldsymbol{s}_{i}(\boldsymbol{z}_{>M})\equiv
a^{i}\boldsymbol{s}_{i}(\boldsymbol{z}_{>M}).$ (17)

As in the second step, Einstein’s convention will be used from now on:
whenever there are repeated indices, the sum $\sum_{1}^{M}$ is implied.

We are mainly be interested in the floating mean periodogram, where

$\displaystyle\boldsymbol{s}_{1}(t)$ $\displaystyle=\sin 2\pi\,t/P$ (18)
$\displaystyle\boldsymbol{s}_{2}(t)$ $\displaystyle=\cos 2\pi\,t/P$
$\displaystyle\boldsymbol{s}_{3}(t)$ $\displaystyle=1$

and $z_{4}=P$ is the period of the signal.

We would like to maximize the log-likelihood ratio between the two hypotheses

$\Delta\chi^{2}=2\log
p(\boldsymbol{x}|\boldsymbol{z},\Sigma)/p(\boldsymbol{x}|\Sigma)=\innerproduct{\boldsymbol{x}}{\boldsymbol{x}}-\innerproduct{\boldsymbol{x}-\boldsymbol{s}}{\boldsymbol{x}-\boldsymbol{s}},$
(19)

where we have introduced a scalar product
$\innerproduct{\boldsymbol{x}}{\boldsymbol{y}}=\boldsymbol{x}^{T}\Sigma^{-1}\boldsymbol{y}$.
Let’s define the metric

$g_{ij}=\innerproduct{\boldsymbol{s}_{i}}{\boldsymbol{s}_{j}}$ (20)

and its inverse $g^{ij}$, such that $g^{ij}g_{jk}=\delta^{i}_{\,k}$. We only
consider problems with a low number of linear parameters, so inverting the
metric is cheap. Let’s define the components of the data along the signal
vectors as

$x_{i}=\innerproduct{\boldsymbol{x}}{\boldsymbol{s}_{i}}.$ (21)

With this notation the score is

$\Delta\chi^{2}=2a^{i}x_{i}-a^{i}a^{j}g_{ij}.$ (22)

Maximization over the linear parameters can be performed analytically:

$\frac{\partial\Delta\chi^{2}}{\partial
a^{i}}\biggr{\rvert}_{\widehat{a}}=2x_{i}-2g_{ij}\widehat{a}^{i}=0,$ (23)

which gives the optimal amplitudes

$\widehat{a}^{i}=g^{ij}x_{j}.$ (24)

The score $\Delta\chi^{2}$ at the optimal parameters is

$\Delta\chi^{2}=2g^{ij}x_{i}x_{j}-g^{ij}g^{kl}x_{j}x_{k}g_{ik}=g^{ij}x_{i}x_{j}.$
(25)

This is the matched filter: the data is matched against the templates
$\boldsymbol{s}_{i}$ and inverse weighted with the noise power, here included
in the definition of the scalar product. Most of the searches to date ignore
the inverse covariance matrix weighting, resulting in a sub-optimal matched
filter analysis.

The strategy for finding the optimal parameters $\boldsymbol{z}$ is to compute
the score from Equation (25) on a grid of non-linear parameters and identify
the highest score.

The main cost of Equation (25) is in computing the scalar products
$\innerproduct{\cdot}{\cdot}$. If noise is stationary with measurements taken
on a regular grid, the computation can be simplified significantly, but we do
not make these assumptions here. The scalar product is computed as
$\innerproduct{\boldsymbol{u}}{\boldsymbol{v}}=\boldsymbol{u}\cdot\widetilde{\boldsymbol{v}}$,
where $\cdot$ is the standard scalar product and $\widetilde{\boldsymbol{v}}$
is the solution of the linear system
$\Sigma\widetilde{\boldsymbol{v}}=\boldsymbol{v}$. Since the covariance matrix
is positive definite we can decompose it using a Cholesky decomposition as
$\Sigma=LL^{T}$, where $L$ is lower triangular. Once the Cholesky
decomposition is computed, solving the linear system is easy: first one solves
$L\boldsymbol{v}^{\prime}=\boldsymbol{v}$ for $\boldsymbol{v}^{\prime}$ and
then $L^{T}\widetilde{\boldsymbol{v}}=\boldsymbol{v}^{\prime}$ for
$\widetilde{\boldsymbol{v}}$.

##  Appendix B Efficient Bayes Factor computation

The difficulty in computing the Bayes Factor (Equation (5)) is in doing the
evidence integrals

$p(X|\mathcal{H}_{i})=\int
p(X|\boldsymbol{z}_{i})p(\boldsymbol{z})d\boldsymbol{z}_{i}.$ (26)

Linear parameters have a Gaussian posterior so can be integrated out of
Equation (26) analytically. They contribute a constant which depends only on
the prior volume of the linear parameters and will therefore be irrelevant if
the Bayes Factor is not taken at a face value, but as a test statistic (Robnik
& Seljak, 2022). What remains is therefore to integrate over the non-linear
parameters at optimal amplitude parameters
$\boldsymbol{a}(\boldsymbol{x},\boldsymbol{z}_{>M})$:

$\displaystyle\widetilde{p}(x|\mathcal{H}_{1})=\int
p(\boldsymbol{x}|\widehat{\boldsymbol{a}}(\boldsymbol{x},\boldsymbol{z}_{>M}),\boldsymbol{z}_{>M})p(\boldsymbol{z}_{>M})d\boldsymbol{z}_{>M}.$
(27)

Because the number of non-linear parameters is low, we can avoid the cost of
MCMC by using the Gaussian quadrature scheme. Our strategy is to find the
parameters which maximize the posterior density and compute the Hessian of the
negative log posterior at those parameters. This is the Lapalace approximation
to the posterior. It is used to set out the Gaussian quadrature scheme and
integrate the posterior (Robnik & Seljak, 2022). In Section 4, the null has
two non-linear parameters and the alternative has one additional non-linear
parameter. We will use quadrature schemes of order 6 and 7 from Stroud &
Secrest (1963). We have compared these computations with pocoMC (Karamanis et
al., 2022a, b), which is a state of the art MCMC method. Both methods agree to
within the MCMC error but quadrature uses 12 log-posterior evaluations in two
dimensions and 27 evaluations in three dimensions, while pocoMC uses $2\times
10^{5}$ in its default setting and still several thousand calls if specialized
to the simple low dimensional targets.

Thus the procedure for computing the Bayes Factor is:

  1. 1.

Maximize the null log-posterior to find the optimal patamers of the null and
the Hessian at the peak.

  2. 2.

Using the Gaussian quadrature scheme, compute the evidence for the null.

  3. 3.

Take the optimal null parameters in the matched filter periodogram from
Section A and find the optimal parameters of the alternative hypothesis.

  4. 4.

Maximize the alternative hypothesis log-posterior over all parameters of the
alternative hypothesis (including the null parameters).

  5. 5.

Using the Gaussian quadrature scheme, compute the evidence for the
alternative.

##  Appendix C Interpolation from the null to the periodic template

Figure 12: Detection count as a function of the Bayes Factor detection
threshold. Left: periodic template and the null-signal template are compared
on real data, as in Figure 8, but with the addition of the $\beta=0.2$ NST.
$\beta=0.2$ template yields a distribution which is practically
indistinguishable from the distribution with the periodic template. Right:
same but with periodic signal injection in the real data, equivalent to Figure
10 in the main text. Even at $\beta=0.2$, NST yields noticably lower scores
than the periodic template, suggesting that detections in the left panel are
false.

In Figure 8 we saw that the NST yields higher scores than the periodic
template. This suggests that there is more non-periodic than periodic features
in the data. As an additional test we here make the NST perturbation from
periodicity less drastic. We introduce a family of null-signal templates,
parametrized by the mixture parameter $\beta$. NSTs in this family are of the
form described in Section 2 and have

$p_{n}=\beta\,p_{n}^{\mathrm{(optimal)}}+(1-\beta),$ (28)

where $p_{n}^{\mathrm{(optimal)}}$ are the values from Section 2. Thus this
family smoothly interpolates between the periodic template ($\beta=0$) and the
NST from Section 2 ($\beta=1$). In Figure 12 we reanalyze the real data
(equivalent of Figure 8 in the main text) and real data with injected signal
(equivalent of Figure 10 in the main text) with NST that has $\beta=0.2$. The
real data analysis shows that the NST gives the same test statistic
distribution as the periodic template. The injected data demonstrates that if
the real signal is present in the data, NST yields lower scores.

## References

  * Baluev (2008) Baluev R. V., 2008, Monthly Notices of the Royal Astronomical Society, 385, 1279
  * Baluev (2013) Baluev R. V., 2013, Monthly Notices of the Royal Astronomical Society, 429, 2052
  * Bayer & Seljak (2020) Bayer A. E., Seljak U., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 009
  * Bayer et al. (2021) Bayer A. E., Seljak U., Robnik J., 2021, Monthly Notices of the Royal Astronomical Society, 508, 1346
  * Charisi et al. (2016) Charisi M., Bartos I., Haiman Z., Price-Whelan A. M., Graham M. J., Bellm E. C., Laher R. R., Márka S., 2016, Monthly Notices of the Royal Astronomical Society, 463, 2145
  * Charisi et al. (2022) Charisi M., Taylor S. R., Runnoe J., Bogdanovic T., Trump J. R., 2022, Monthly Notices of the Royal Astronomical Society, 510, 5929
  * Chen et al. (2020) Chen Y.-C., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 2245
  * Chen et al. (2024) Chen Y.-J., et al., 2024, Monthly Notices of the Royal Astronomical Society, 527, 12154
  * D’Orazio & Charisi (2023) D’Orazio D. J., Charisi M., 2023, Observational Signatures of Supermassive Black Hole Binaries, doi:10.48550/arXiv.2310.16896, https://ui.adsabs.harvard.edu/abs/2023arXiv231016896D
  * D’Orazio et al. (2015) D’Orazio D. J., Haiman Z., Schiminovich D., 2015, Nature, 525, 351
  * Duffell et al. (2020) Duffell P. C., D’Orazio D., Derdzinski A., Haiman Z., MacFadyen A., Rosen A. L., Zrake J., 2020, The Astrophysical Journal, 901, 25
  * D’Orazio & Stefano (2020) D’Orazio D. J., Stefano R. D., 2020, Monthly Notices of the Royal Astronomical Society, 491, 1506
  * Fowlie (2023) Fowlie A., 2023, Communications in Statistics - Theory and Methods, 52, 5379
  * Graham et al. (2015) Graham M. J., et al., 2015, Monthly Notices of the Royal Astronomical Society, 453, 1562
  * Haiman et al. (2009) Haiman Z., Kocsis B., Menou K., 2009, The Astrophysical Journal, 700, 1952
  * Haiman et al. (2023) Haiman Z., et al., 2023, Massive Black Hole Binaries as LISA Precursors in the Roman High Latitude Time Domain Survey, doi:10.48550/arXiv.2306.14990, http://arxiv.org/abs/2306.14990
  * Ivezic et al. (2019a) Ivezic Z., et al., 2019a, The Astrophysical Journal, 873, 111
  * Ivezic et al. (2019b) Ivezic Z., et al., 2019b, The Astrophysical Journal, 873, 111
  * Jeffreys (1998) Jeffreys H., 1998, The Theory of Probability. OUP Oxford
  * Jon M. Jenkins et al. (2017) Jon M. Jenkins Tenenbaum Seader Burke McCauliff Smith Twicken Chandrasekaran 2017, Kepler Data Processing Handbook: Transiting Planet Search. https://ui.adsabs.harvard.edu/abs/2017ksci.rept....9J/abstract
  * Karamanis et al. (2022a) Karamanis M., Beutler F., Peacock J. A., Nabergoj D., Seljak U., 2022a, Astrophysics Source Code Library, p. ascl:2207.018
  * Karamanis et al. (2022b) Karamanis M., Beutler F., Peacock J. A., Nabergoj D., Seljak U., 2022b, Monthly Notices of the Royal Astronomical Society, 516, 1644
  * Kelly et al. (2009) Kelly B. C., Bechtold J., Siemiginowska A., 2009, The Astrophysical Journal, 698, 895
  * Kozłowski et al. (2010) Kozłowski S., et al., 2010, The Astrophysical Journal, 708, 927
  * LIGO Scientific Collaboration et al. (2004) LIGO Scientific Collaboration et al., 2004, Physical Review D, 69, 122001
  * Law et al. (2009) Law N. M., et al., 2009, Publications of the Astronomical Society of the Pacific, 121, 1395
  * Liu et al. (2018) Liu T., Gezari S., Miller M. C., 2018, The Astrophysical Journal, 859, L12
  * Liu et al. (2019) Liu T., et al., 2019, The Astrophysical Journal, 884, 36
  * Lomb (1976) Lomb N. R., 1976, Astrophysics and Space Science, 39, 447
  * MacLeod et al. (2010) MacLeod C. L., et al., 2010, The Astrophysical Journal, 721, 1014
  * Peters (1964) Peters P. C., 1964, Physical Review, 136, B1224
  * Robnik & Seljak (2021) Robnik J., Seljak U., 2021, Monthly Notices of the Royal Astronomical Society, 504, 5829
  * Robnik & Seljak (2022) Robnik J., Seljak U., 2022, Entropy, 24, 1328
  * Scargle (1982) Scargle J. D., 1982, The Astrophysical Journal, 263, 835
  * Stroud & Secrest (1963) Stroud A. H., Secrest D., 1963, Mathematics of Computation, 17, 105
  * VanderPlas (2018) VanderPlas J. T., 2018, The Astrophysical Journal Supplement Series, 236, 16
  * Vaughan et al. (2016) Vaughan S., Uttley P., Markowitz A. G., Huppenkothen D., Middleton M. J., Alston W. N., Scargle J. D., Farr W. M., 2016, Monthly Notices of the Royal Astronomical Society, 461, 3145
  * Weakliem (1999) Weakliem D. L., 1999, Sociological Methods & Research, 27, 359
  * Westernacher-Schneider et al. (2022) Westernacher-Schneider J. R., Zrake J., MacFadyen A., Haiman Z., 2022, Physical Review D, 106, 103010
  * Witt et al. (2022) Witt C. A., Charisi M., Taylor S. R., Burke-Spolaor S., 2022, The Astrophysical Journal, 936, 89
  * Xin & Haiman (2021) Xin C., Haiman Z., 2021, Monthly Notices of the Royal Astronomical Society, 506, 2408
  * Zhang (2017) Zhang J., 2017, Australian & New Zealand Journal of Statistics, 59, 43
  * Zhu & Thrane (2020) Zhu X.-J., Thrane E., 2020, The Astrophysical Journal, 900, 117