# Pitfalls of periodograms: The non-stationarity bias in the analysis of
quasi-periodic oscillations

Moritz Hübner<EMAIL_ADDRESS>School of Physics and Astronomy,
Monash University, Clayton, VIC 3800, Australia  OzGrav: The ARC Centre of
Excellence for Gravitational Wave Discovery, Clayton, VIC 3800, Australia
Daniela Huppenkothen  SRON Netherlands Institute for Space Research  Paul D.
Lasky  School of Physics and Astronomy, Monash University, Clayton, VIC 3800,
Australia  OzGrav: The ARC Centre of Excellence for Gravitational Wave
Discovery, Clayton, VIC 3800, Australia  Andrew R. Inglis  Physics Department,
The Catholic University of America, Washington, DC, 20064, USA  Solar Physics
Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, 20771, USA

###### Abstract

Quasi-periodic oscillations (QPOs) are an important key to understand the
dynamic behavior of astrophysical objects during transient events like gamma-
ray bursts, solar flares, and magnetar flares. Searches for QPOs often use the
periodogram of the time series and perform spectral density estimation using a
Whittle likelihood function. However, the Whittle likelihood is only valid if
the time series is stationary since the frequency bins are otherwise not
statistically independent. We show that if time series are non-stationary, the
significance of QPOs can be highly overestimated and estimates of the central
frequencies and QPO widths can be overconstrained. The effect occurs if the
QPO is only present for a fraction of the time series and the noise level is
varying throughout the time series. This can occur for example if background
noise from before or after the transient is included in the time series or if
the low-frequency noise profile varies strongly over the time series. We
confirm the presence of this bias in previously reported results from solar
flare data and show that significance can be highly overstated. Finally, we
provide some suggestions that help identify if an analysis is affected by this
bias.

##  1 Introduction

Quasi-periodic oscillations (QPOs)111In solar physics the term “quasi-periodic
pulsation” (QPP) is preferred. We apply the term QPO generically throughout
for simplicity. are a common astrophysical phenomenon that are regularly
observed across a variety of sources. While there is ample discussion about
how they emerge in their respective sources, it is worth re-examining existing
techniques for their detection and characterization. Observations of QPOs are
scant for some objects such as in magnetar flares (Strohmayer & Watts, 2006a,
2005; Israel et al., 2005; Watts & Strohmayer, 2006; Huppenkothen et al.,
2012, 2014b, 2014a; Miller et al., 2019), or contentious in others like gamma-
ray bursts (GRBs) (Cenko et al., 2010; De Luca et al., 2010; Iwakiri et al.,
2010; Morris et al., 2010; Tarnopolski & Marchenko, 2021), and false
detections may lead to unjustified theoretical inferences. Moreover, even if
there are ample detections of QPOs like in solar flares (Nakariakov &
Melnikov, 2009; Van Doorsselaere et al., 2016; Zimovets et al., 2021) or X-ray
binaries (Ingram & Motta, 2019), making accurate inferences about their
parameters and significance is essential for understanding the mechanisms that
produce them.

There is an array of methods being applied to time series data to find QPOs.
The most common tests involve periodograms and are either based on outlier
analyses or Bayesian tests (Vaughan, 2005, 2009). Other methods include
wavelet analyses (Foster, 1996; Torrence & Compo, 1998; Lachowicz & Done,
2010) and those based on Gaussian processes (Rasmussen & Williams, 2006;
Foreman-Mackey et al., 2017; Covino et al., 2020; Zhu & Thrane, 2020).

In this paper, we explore the effects that the non-stationary nature of a
transient lightcurve, and the potential non-stationarity of a QPO within it,
have on detection significance and characterization when methods are used that
assume the underlying processes to be stationary. We show that we are likely
to overestimate the significance of the quasi-periodic signal under these
conditions. This bias arises because neighboring bins in periodograms are only
statistically independent for stationary time series. We thus refer to this
effect as the “non-stationarity bias”, which to the best of our knowledge has
not been reported in the astrophysical literature to date. We show that merely
adding photon counting noise before or after a transient leads one to greatly
overestimate the significance of a QPO and overconstrain the QPO parameters.
More critically, this effect also occurs if only the transient is selected for
analysis but the QPO is only present for a fraction of it.

It is our intuition that astrophysicists would rather be conservative and
include more of the time series that does not include the QPO, rather than
“cherry-picking” the part of the lightcurve that appears most likely to
contain it. The erroneous reasoning behind this might be the belief that
adding more noise to the signal should not increase the significance of the
QPO.

We structure this paper as follows. In Sec. 2 we lay out the Bayesian methods
and models that we use throughout. Thereafter, in Sec. 3, we explain how the
non-stationarity bias arises both on a conceptual level and with detailed
mathematical arguments. We show in Sec. 4 based on simulated data that the
non-stationarity bias exists empirically. In Sec. 5 we show how the non-
stationarity bias affects the analysis of solar flare lightcurves. We conclude
in Sec. 6 and provide some outlook on how alternative methods can potentially
better handle non-stationary data sets.

##  2 Methods

In the following, we recapitulate the Fourier-based methods commonly used to
analyze QPOs which are relevant to this study. A general comprehensive
overview on the topic of spectral density estimation for time series analysis
is given in many popular textbooks such as chapter 7 in Chatfield (1975), as
well as in van der Klis (1989); Barret & Vaughan (2012) in the context of
astrophysical lightcurves.

###  2.1 Whittle likelihood

A periodogram is an estimate of the power spectral density (PSD) of the signal
based on a given time series $x(t)$. We calculate the periodogram $I(f)$ as
the absolute square of the discrete Fourier transform of the signal

$I(f)=\left|\sum_{i=1}^{N}x(t_{i})\exp(-i2\pi ft_{i})\right|^{2}\,.$ (1)

Periodograms are established as a standard method in the search for QPOs in
solar flares (Inglis et al., 2015, 2016; Hayes et al., 2020) and QPOs in
astrophysical transients elsewhere (e.g. Huppenkothen et al., 2012, 2014b,
2014a, 2017; Miller et al., 2019). Since individual frequency bins
$I(f_{j})=I_{j}$ are calculated by taking the sum of the square of the normal-
distributed real and imaginary parts of the Fourier series, it follows that
they are $\chi^{2}_{2}$-distributed around the PSD $S(f_{j})=S_{j}$ (Whittle,
1951; Chatfield, 1975). This specific case of a $\chi^{2}_{2}$-distribution is
identical to an exponential distribution around $S_{j}$

$p(I_{j}|S_{j})=\frac{1}{S_{j}}\exp(-I_{j}/S_{j})\,,$ (2)

where $p(I_{j}|S_{j})$ is the conditional probability to observe the power
$I_{j}$ given an underlying PSD $S_{j}$. This relation is generically true for
any individual frequency bin in any periodogram. Assuming that all bins are
statistically independent from another, we obtain the Whittle likelihood
function by taking the product over all $N/2$ frequency bins of a periodogram
corresponding to a time series with $N$ points

$L(I|S)=\prod_{j=1}^{N/2}\frac{1}{S_{j}}\exp(-I_{j}/S_{j})\,.$ (3)

We emphasize here that this is only true in the stationary limit, and this
will lead to biased estimates of the PSD in general, as we show in Sec. 3.
Nevertheless, Eq. 3 is the standard likelihood that is used for spectral
density estimation.

###  2.2 Models of power spectra

Many astrophysical transients show excess power at low frequencies, and it is
often assumed that this can be modeled using a red-noise process (see e.g.
Huppenkothen et al. 2012; Inglis et al. 2015, 2016; Miller et al. 2019; Ingram
& Motta 2019; Broomhall et al. 2019). One basic noise model that is commonly
used is a combination of a red-noise power law with amplitude $A$ and spectral
index $\alpha$, and a white noise amplitude $C$ (Inglis et al., 2015, 2016),
i.e.

$\displaystyle S^{\mathrm{RW}}(f)$
$\displaystyle=S^{\mathrm{R}}(f)+S^{\mathrm{W}}$ (4)
$\displaystyle=Af^{-\alpha}+C\,,$

where we use the $\mathrm{R}$ and $\mathrm{W}$ superscripts as short-hand for
“red noise” and “white noise” respectively. This spectral shape usually
emerges as a combination from the overall structure of the transient, and
additional variability on smaller time scales. The shape is thus not due to a
stationary process but rather due to a combination of unknown deterministic
and non-stationary stochastic processes adding up to mimic a red noise
spectrum (Huppenkothen et al., 2012).

A common way to model a QPO is to add a Gaussian or Lorentzian enhancement to
$S^{\mathrm{RW}}$, e.g.

$\displaystyle S^{\mathrm{RWQ}}(f)$
$\displaystyle=S^{\mathrm{RW}}(f)+S^{\mathrm{Q}}(f)$ (5)
$\displaystyle=S^{\mathrm{RW}}(f)+\frac{B}{\pi\sigma}\frac{\sigma^{2}}{(f-f_{0})^{2}+\sigma^{2}}\,,$

where $B$ is the QPO amplitude, $f_{0}$ is its mean frequency, and $\sigma$ is
the half-width half-maximum scale parameter. Explicitly modeling the red noise
component is important as the QPO is likely to overlap with the red noise
dominated part of the PSD. Not considering the red noise dominated part of the
PSD may thus lead to false positives.

There are extensions to the noise model, such as the broken or bent power law
(B superscript), that can fit more structured red noise and are often a better
fit to the data

$S^{\mathrm{BW}}(f)=Af^{-\alpha_{1}}\left(1+\left(\frac{f}{\delta}\right)^{\alpha_{2}-\alpha_{1}}\right)+C\,,$
(6)

where $\alpha_{1,2}$ are the power law indices present before and after the
break frequency $\delta$ where the power law changes, and we also enforce
$\alpha_{2}<\alpha_{1}$ to avoid degeneracies. There are also various
different formulations of bent or broken power laws with various degrees of
smoothness.

###  2.3 Model selection and parameter estimation

There are several ways to assert the significance of a QPO when using spectral
density estimation. A widely used frequentist way to detect QPOs is to use
outlier statistics. As a first step, one fits the PSD using a model that does
not contain a QPO, e.g. Eq. 4 or Eq. 6. Thereafter, every frequency bin is
checked if its amplitude exceeds a set probability threshold based on the
$\chi^{2}_{2}$-distribution in Eq. 2. One also has to account for the number
of trials, i.e. the number of frequency bins tested by applying a Bonferroni-
correction (Bonferroni, 1936). Additionally, it is possible to rebin the
periodogram into a smaller number of bins. For QPOs with widths larger than
the width of a single frequency bin, rebinning may be useful because it
increases the significance of the QPO relative to the noise continuum.

Alternative to the outlier statistics, one can cast this as a model Bayesian
selection problem where we find the preferred model to fit our data, in our
case Eqs. 4, 5, and 6. There are multiple ways to quantify model preference in
Bayesian statistics. One approach that has been used in solar physics (Ireland
et al., 2015; Inglis et al., 2015, 2016) is to calculate the Bayesian
Information Criterion ($BIC$)

$BIC=k\ln(n)-2\ln(L_{\max})$ (7)

where $k$ is the number of free parameters (e.g. three for $S^{\mathrm{RW}}$,
six for $S^{\mathrm{RWQ}}$, and five for $S^{\mathrm{BW}}$), $n=N/2$ is the
number of data points, and $L_{\max}$ is the maximum likelihood value. This
method is relatively cheap computationally as the calculation of the maximum
likelihood can be obtained with relatively few likelihood evaluations using a
maximization algorithm. A lower $BIC$ indicates a better fit to the data, thus
the difference in $BIC$ for $S^{\mathrm{RW}}$ and $S^{\mathrm{RWQ}}$ is a
measure of statistical significance of the QPO,

$\displaystyle\Delta BIC$
$\displaystyle=BIC(S^{\mathrm{RWQ}})-BIC(S^{\mathrm{RW}})$ (8)
$\displaystyle=(k_{\mathrm{RWQ}}-k_{\mathrm{RW}})\ln(N/2)-2\ln\frac{L_{\max}(I|S^{\mathrm{RWQ}})}{L_{\max}(I|S^{\mathrm{RW}})}\,,$

where $\Delta BIC<0$ would indicate that $S^{\mathrm{RWQ}}$ is preferred and
vice versa. We can also perform model selection via calculation of Bayes
factors $BF$. To understand Bayes factors, we start with Bayes’ theorem

$p(\theta|d,S)=\frac{\pi(\theta|S)L(d|\theta,S)}{Z(d|S)}\,,$ (9)

where $\theta$ are the model parameters, $d$ are the data, i.e. the
periodogram in our case, $p$ is the posterior probability of the parameters,
$\pi$ is the prior probability of the parameters, $L$ is the likelihood of the
data given the parameters, and $Z$ is the evidence, or fully marginalized
likelihood. All these probabilities are conditioned on a model we want to
evaluate, e.g. a PSD $S$. The Bayes factor comparing two models is the ratio
of their evidences. For example, the Bayes factor comparing $S^{\mathrm{RWQ}}$
and $S^{\mathrm{RW}}$ is

$BF=\frac{Z(d|S^{\mathrm{RWQ}})}{Z(d|S^{\mathrm{RW}})}\,.$ (10)

It thereby measures the relative odds of the underlying data to have been
produced by either model, though it does not measure if the model itself is a
good fit to the data, similar to the $BIC$. The evidence is calculated by
rearranging and integrating Bayes’ theorem

$Z(d|S)=\int\pi(\theta)L(d|\theta,S)d\theta\,.$ (11)

Evaluating this integral is much more computationally challenging than
calculating the $BIC$, but can be achieved thanks to improvements in
algorithms such as nested sampling (Skilling, 2006), and accessible software
implementations such as Bilby (Speagle, 2019; Ashton et al., 2019; Romero-Shaw
et al., 2020). The Bayes factor obtained via evidence calculation is seen as
the superior standard for model selection in Bayesian statistics since it
involves prior beliefs about the distribution of parameters (Weakliem, 1999).
The penalty factor $k\ln(n)$ for the $BIC$ also increases with the size of the
data set and can thus favor overly simplistic models (Weakliem, 1999; Gelman
et al., 2013). Since the $BIC$ relies on a point estimate of the likelihood,
it is also a less reliable measure if the likelihood has multiple modes.
Additionally, nested sampling provides posterior distributions on the model
parameters which are interesting in their own right.

Bayes factors and $BIC$s work best if the underlying models are truly
discrete, but $S^{\mathrm{RW}}$ is a special case of both $S^{\mathrm{RWQ}}$
and $S^{\mathrm{BW}}$. There are reasonable criticisms on the use of Bayes
factors or the $BIC$ for cases like this as $S^{\mathrm{RW}}$ is a subset of
both $S^{\mathrm{RWQ}}$ and $S^{\mathrm{BW}}$, and thus not a truly discrete
model (Gelman et al., 2013). In this case, the Bayes factor depends on prior
choices of the parameters specific to $S^{\mathrm{RWQ}}$ and
$S^{\mathrm{BW}}$. Alternatively, we could use tests that measure if the
posterior of the QPO amplitude $B$ is inconsistent with being zero and
provides a better fit to the data. Nevertheless, we use Bayes factors as prior
choices since they will yield a constant difference in our Bayes factors and
will not distort the overall trend when we investigate the non-stationarity
bias. We also consider the effects on the $\Delta BIC$s where it is
instructive.

For all the Bayes factor calculations and posterior samples in this
publication we use dynesty via the Bilby interface (Ashton et al., 2019;
Speagle, 2019; Romero-Shaw et al., 2020).

##  3 Pitfalls of Periodograms

Figure 1: Effects of extending a time series with zeros. In the top panel, we
show a stationary time series between $-10\text{\,}\mathrm{s}$ and
$10\text{\,}\mathrm{s}$ that we multiply with a Hann window to ensure a smooth
turn on from zero. We extended this time series with zeros on either end up to
a total duration of $200\text{\,}\mathrm{s}$. The upper x-axis shows how the
extension factor $x$ maps onto the time axis. Taking only the stationary data
from $-10\text{\,}\mathrm{s}$ to $10\text{\,}\mathrm{s}$ corresponds to $x=1$,
i.e. no extension has been applied. The extension factor in this case implies
that we select data from $-10x$ $\mathrm{s}$ to $10x$ $\mathrm{s}$. In the
bottom panel, we show the periodogram given different extension factors, with
the blue curve representing the stationary data ($x=1$). The other colors
(orange, green, red) show the effect of extending the time series with zeros
on either end. In these curves, neighboring bins are not statistically
independent, which manifests itself in the emerging smooth structure of the
periodogram. The bins are also decreasing inversely proportional to $x$ since
we do not apply a normalization.

Assume a discrete time series of length $T$ with $N$ data points containing a
QPO and noise. The underlying stochastic process is characterized in the
frequency domain as a PSD $S(f)$, which can be estimated by calculating a
periodogram $I(f)$. We recall that the value in any specific of the $N/2$
frequency bins in the periodogram is proportional to the absolute square of
the Fourier amplitude.

In the top panel of Fig. 1 we extend a stationary time series by appending
zeros on either end until we reach a length $xT$ with $xN$ data points, where
we call $x$ the extension factor. As we display in the bottom panel, this
increases the number of frequency bins in the periodogram to $xN/2$, while
decreasing the periodogram to $I(f)/x$ if we do not apply a normalization. The
non-stationarity bias occurs here because an increased number of frequency
bins boosts the calculated significance of a QPO in terms of the $\Delta BIC$
or $\ln BF$ between, e.g. $S^{\mathrm{RWQ}}$ and $S^{\mathrm{RW}}$.

The example displayed in Fig. 1 is not realistic in practice but serves to
illustrate the effect. A more realistic scenario occurs if the stationary time
series is instead extended with low amplitude white noise and the QPO is
within the red noise dominated part of the PSD. This raises the white noise
level relative to red noise and the QPO. Yet it still increases the
significance of the QPO, which is determined by its amplitude relative to the
dominant red noise level. Such a scenario might occur if we overselect data
surrounding a transient

Another scenario is that the overall time series is non-stationary, and
instead, the QPO and the overall noise processes vary, or switch on and off,
during a selected time window. For example, the QPOs reported in the 2004
SGR1806-20 hyperflare are directly associated with specific rotational phases
of the magnetar (Strohmayer & Watts, 2005; Israel et al., 2005; Strohmayer &
Watts, 2006b; Watts & Strohmayer, 2006; Huppenkothen et al., 2014a; Miller et
al., 2019). Parts of the time series might only contribute low levels of noise
and thereby increase the number of bins without suppressing the QPO in the
periodogram. These parts, therefore, increase the significance of the QPO
without having delivered any actual information about it.

We now turn towards a more mathematical explanation of the non-stationarity
bias. Let us consider how we can calculate the SNR of a QPO when red and white
noise are present222The argument here holds for arbitrary noise.. Given a
single frequency bin $f_{j}$, this is

$\rho(f_{j})=\frac{S^{\mathrm{Q}}(f_{j})}{S^{\mathrm{RW}}(f_{j})}\,,$ (12)

where we denote the SNR as $\rho$. Assuming all frequency bins in a
periodogram are statistically independent, the SNR adds in quadrature

$\rho_{\mathrm{tot}}=\sqrt{\sum_{j=1}^{N/2}\rho(f_{j})^{2}}\,,$ (13)

where we have summed over the $N/2$ bins in the periodogram and
$\rho_{\mathrm{tot}}$ denotes the total SNR of the QPO.

###  3.1 QPOs in stationary noise

As the most simple scenario, we consider a stationary time series with $N$
elements made up of a QPO and some arbitrary noise. We denote the SNR of this
time series as $\rho^{N}_{\mathrm{tot}}$. Suppose we look at an extended
version of this time series with $xN$ elements, which we create with the same
QPO and noise process. We intuitively expect the SNR to increase
proportionally to $\sqrt{x}$ and indeed, we calculate

$\displaystyle\begin{split}\rho^{xN}_{\mathrm{tot}}&=\sqrt{\sum_{j=1}^{xN/2}\rho^{x}(f_{j})^{2}}\\\
&=\sqrt{x\sum_{j=1}^{N/2}\rho(f_{j})^{2}}\\\
&=\sqrt{x}\rho^{N}_{\mathrm{tot}}\,,\end{split}$ (14)

where we write the SNR at a given frequency in the extended time series as
$\rho^{x}(f_{j})$, which is the same as in the original time series. We also
use the $N$ and $xN$ superscripts for the SNRs to indicate whether we are
looking at the original or extended time series. Note that we have implicitly
changed how the $f_{j}$ are indexed between the steps. Additionally, we assume
in the first step that the PSD is only slowly varying between frequency bins
and we can therefore approximate the newly added frequency bins by their
closest neighbors from the original PSD.

Alternatively, let us consider what happens when we consider a non-stationary
extended time series in which the noise is present over $xN$ elements but the
QPO is only present for $N$ elements. In that case, $S^{\mathrm{RW}}$ is the
same as in the original time series whereas $S^{\mathrm{Q}}$ gets reduced to
$S^{\mathrm{Q}}/x$ as we are effectively “diluting” the QPO, similarly to how
the periodograms decrease in amplitude with increasing $x$ in the lower panel
of Fig. 1. We calculate Eq. 13 for the extended time series again

$\displaystyle\begin{split}\rho^{xN}_{\mathrm{tot}}&=\sqrt{\sum_{j=1}^{xN/2}\rho^{x}(f_{j})^{2}}\\\
&=\sqrt{x\sum_{j=1}^{N/2}\rho^{x}(f_{j})^{2}}\\\
&=\sqrt{x\sum_{j=1}^{N/2}\rho(f_{j})^{2}/x^{2}}\\\
&=\frac{1}{\sqrt{x}}\rho^{N}_{\mathrm{tot}}\,,\end{split}$ (15)

where we used the first same step as in Eq. 14 and used the fact that the QPO
is suppressed by the extension in the second step. These first two scenarios
are fairly intuitive, analyzing more of the same stationary process should
naturally boost the significance, whereas only adding more noise will dilute
the QPO and reduce the significance.

Let us now consider what happens if we extend the time series with zeros on
either end, which is illustrated in Fig. 1. This scenario is not truly
physical, as we always expect some background noise. However, it is an
instructive example to understand the effects of non-stationarities. Extending
with zeros implies both $S^{\mathrm{RWQ}}\rightarrow S^{\mathrm{RWQ}}/x$ and
$S^{\mathrm{RW}}\rightarrow S^{\mathrm{RW}}/x$, so $\rho^{x}(f)=\rho(f)$ and
thus following the same steps as before

$\rho^{xN}_{\mathrm{tot}}=\sqrt{x}\rho^{N}_{\mathrm{tot}}\,,$ (16)

which is identical to Eq. 14. This is a curious result as extending a time
series by zeros clearly should not increase the SNR nor the significance of a
signal. The resolution to this seeming paradox is that it is invalid to assume
that SNR adds in quadrature as we did in Eq. 14. By extending the time series
with zeros, we added bins to the periodogram that are not statistically
independent.

A further consequence is that the product in the Whittle likelihood becomes
invalid. Due to the extended segment length both the PSD and the periodogram
need to be divided by a factor of $x$. Therefore, we can rewrite the
likelihood as

$\displaystyle L^{xN}(I|S(f|\theta))$
$\displaystyle=\prod_{j=1}^{xN/2}\frac{x}{S_{j}(\theta)}\exp(-I_{j}/S_{j}(\theta))$
(17) $\displaystyle\approx
x^{xN/2}\left(\prod_{j=1}^{N/2}\frac{1}{S_{j}(\theta)}\exp(-I_{j}/S_{j}(\theta))\right)^{x}\,,$

where $\theta$ is the set of parameters in the PSD and we use the $N$ and $xN$
superscipts on the likelihoods to indicate if we are looking at the original
or extended time series. Again, we assume that the PSD is slowly varying with
frequency and thus statistically non-independent bins are reasonably
approximated by the closest frequency bins obtained from the stationary time
series. The factor $x^{xN/2}$ that we have introduced does not matter since it
is a constant. Hence we can write

$\ln L^{xN}(I|S(\theta))\propto x\ln L^{N}(I|S(\theta))\,.$ (18)

This means $\ln L^{xN}(I|S(\theta))\approx x\ln L^{N}(I|S(\theta))$, i.e. we
have steepened the log likelihood function by a factor of $x$

$\frac{\partial\ln L^{xN}(I|S(\theta))}{\partial\theta}\approx
x\frac{\partial\ln L^{N}(I|S(\theta))}{\partial\theta}\,.$ (19)

This steepening means that inferred posterior distributions are generally
tighter which leads to overconstrained parameter estimates. We calculate how
this changes the Bayesian information criterion

$\displaystyle\begin{split}\Delta
BIC^{xN}&=(k_{\mathrm{RWQ}}-k_{\mathrm{RW}})\ln(xN/2)-\\\
&2x\ln\frac{L^{N}_{\max}(I|S^{\mathrm{Q}})}{L^{N}_{\max}(I|S^{\mathrm{RW}})}\\\
&=x\Delta BIC+\\\ &\,(k_{\mathrm{RWQ}}-k_{\mathrm{RW}})(\ln
x-(x-1)\ln(N/2))\,,\end{split}$ (20)

where $k_{\mathrm{RWQ}}=6$ and $k_{\mathrm{RW}}=3$ are the number of free
parameters for the respective models, and $xN$ is used as a superscript. Thus,
$\Delta BIC$ is approximately proportional to $x$ or $N$ since the impact of
the $\ln x$ and constant terms are minimal. The impact this has on the Bayes
factor is not straight-forward to calculate in the general case. Since
likelihood ratios and $\Delta BIC$s can be understood as a related quantity of
the $\ln BF$, we expect a similar approximately linear behavior, which we find
empirically in Sec. 4.

###  3.2 QPOs in non-stationary red noise

Assume a stationary time series with red noise, white noise, and a QPO that
sits in the red-noise dominated part of the PSD. If we extend this time series
with zeros, the argument in Sec. 3.1 holds and the SNR calculated with Eq. 13
grows proportional to $\sqrt{x}$. Extending with white noise instead of zeros
will also lead to the SNR growing with $\sqrt{x}$ as long as the QPO remains
in the red noise dominated part of the PSD. To show this we define
$f_{\mathrm{break}}$ as the breaking frequency where white noise and red noise
are of equal magnitude

$f_{\mathrm{break}}=\left(\frac{C}{A}\right)^{-1/\alpha}\,,$ (21)

which means that for $f\ll f_{\mathrm{break}}$ red noise dominates whereas
white noise dominates for $f\gg f_{\mathrm{break}}$. If our QPO has
$f_{0}<f_{\mathrm{break}}$ it mainly competes with red noise in terms of SNR,
thus a rising level of white noise will not meaningfully reduce the SNR until
$f_{0}\approx f_{\mathrm{break}}$. We can write this as

$\rho_{\mathrm{tot}}\approx\rho_{\mathrm{tot}}(f<f_{\mathrm{break}})\,,$ (22)

where $\rho_{\mathrm{tot}}(f<f_{\mathrm{break}})$ is the total SNR just based
on frequencies less than the breaking frequency. If we extend the time series
with white noise or zeros by a factor of $x$, we suppress the red noise term
by $1/x$ and shift the breaking frequency

$f_{\mathrm{break}}^{x}=\left(\frac{Cx}{A}\right)^{-1/\alpha}=x^{-1/\alpha}f_{\mathrm{break}}\,,$
(23)

where the superscript $x$ indicates that we are looking at the extended time
series. We can equate Eq. 23 to $f_{0}$ and re-arrange to calculate the
extension factor $x_{\mathrm{break}}$ for which the QPO will be in the white-
noise dominated part of the periodogram

$x_{\mathrm{break}}=\left(\frac{f_{\mathrm{break}}}{f_{0}}\right)^{\alpha}=\frac{A}{C}f_{0}^{-\alpha}\,.$
(24)

Therefore, we are prone to artificially increase the SNR by extending our time
series with white noise up to an extension factor of $x_{\mathrm{break}}$.
Conversely, if our QPO mainly competes with white noise
$f_{0}>f_{\mathrm{break}}$, we use the argument in Sec. 3.1 to show that we
suppress the SNR with increasing $x$

$\displaystyle\begin{split}\rho_{\mathrm{tot},x}&=\sqrt{\sum_{j=1}^{xN/2}\left(\frac{S^{\mathrm{Q}}(f_{j})/x}{S^{\mathrm{W}}}\right)^{2}}\\\
&\approx\frac{1}{\sqrt{x}}\sqrt{\sum_{j=1}^{N/2}\left(\frac{S^{\mathrm{Q}}(f_{j})}{S^{\mathrm{W}}}\right)^{2}}\\\
&\approx\frac{\rho_{\mathrm{tot}}}{\sqrt{x}}\,,\end{split}$ (25)

where we have used the fact that the PSD is slowly varying with frequency in
the first approximation, and that the frequencies where the QPO exceeds the
noise continuum are all white-noise dominated in the second approximation.

The most obvious case where this could become a practical problem is when we
have to decide where to start and end the segment. Not being aware of this
bias may lead us to overselect the data by $\mathcal{O}(10\%)$ and infer a
somewhat higher SNR and significance, which by itself is not a big issue. But
the issue becomes much worse in practice if the QPO itself is non-stationary
and only appears for a part of the transient, which implies that selecting the
entire transient may already be overextending the segment in which the QPO is
present.

In principle, extending the time series with red noise instead of white noise
should not cause the same issue as we would expect the same suppression by
$1/\sqrt{x}$ that we have for white noise in Sec. 3.1. However, in real
transients, the low-frequency noise continuum arises not due to a stationary
noise process but due to a combination of deterministic and non-stationary
stochastic processes. Different parts of the time series might add up roughly
to a single power law or a broken power law, but smaller segments within the
transient may have vastly different shapes. For example, the rising and
falling edge of the transient may create different power laws, we may have
segments that are fairly flat and thus mostly add white noise, and so on.
Additionally, deterministic aspects of the time series, if not properly
subtracted, are not appropriately modeled with a Whittle likelihood as we show
in Sec. 4.3. The deterministic parts of the lightcurve can not be safely
removed by filtering or smoothing methods as they are prone to create an
artificial structure that looks like oscillatory behavior (Auchère et al.,
2016). Clearly, it is not possible to write down a rigorous treatment of this
rather large class of possible constellations in which non-stationarity bias
could occur.

We note here that using the outlier analysis is not prone to false detections
due to the non-stationarity bias since we are testing neighboring bins
individually instead of combining SNR from neighboring bins. The opposite is
the case, since we add more bins by having extended the time series, it is
harder to exceed the Bonferroni-corrected significance level. However, in
practice periodograms are sometimes rebinned to add up neighboring bins (e.g.
Huppenkothen et al., 2014b). The motivation behind this is to combine the SNR
from neighboring bins, making it easier to detect QPOs. Thus, rebinned outlier
analyses are prone to the same non-stationarity bias that affects the Bayesian
methods as well.

##  4 Simulated data

To demonstrate the impact of non-stationarities empirically, we analyze
simulated data using a variety of setups. For all studies on simulated data,
we assume white noise levels to be constant along with the time series. We
discuss the impact of the Poissonian nature of the photon counting process in
Appendix B.

###  4.1 Setup

To produce the time domain data, we use the algorithm by Timmer & König
(1995), as it is implemented in the Bilby software package. Concretely, this
works by creating a white noise frequency series by randomly drawing both
amplitudes and phases from a normal distribution, coloring the noise by
multiplying it with the square root of the PSD and then applying an inverse
Fourier transform to obtain a stationary time series. In general, we create a
continuous, non-stationary time series by using a combination of addition,
concatenation, and convolution with window functions. For example, we start by
creating a stationary time series using the QPO model and red noise
$S^{\mathrm{RQ}}=S^{\mathrm{R}}+S^{\mathrm{Q}}$. We apply a Hann window to
ensure a smooth turn on from zero. This helps us to avoid Fourier artifacts
such as side-lobes in the periodogram which can appear if discontinuities are
present. Next, we create a much longer time series containing just white noise
using the Timmer & König (1995) method, and add the $S^{\mathrm{RQ}}$ time
series in the center. We are also interested in what happens when we extend
the time series with zeros instead of white noise. To do this we mask the
parts containing just white noise and set those values to zero. In all cases,
we center the part of the time series in which the QPO is present. This is to
prevent the Hann window, which is applied when calculating the periodogram,
from diminishing the amplitude of the red noise and QPO components. Note that
this is separate from the Hann window we use to ensure the smooth turn-on from
zero.

For all simulated data in this section, we assume a constant white noise level
throughout the time series and place the QPO in the center of the time series
between $-10\text{\,}\mathrm{s}$ and $10\text{\,}\mathrm{s}$. We use a
sampling frequency of $40\text{\,}\mathrm{Hz}$ and extend the time series up
to $200\text{\,}\mathrm{s}$ ($x=10$) or $400\text{\,}\mathrm{s}$ ($x=20$)
depending on the scenario. The specific parameters for each simulation are
listed in the tables in App. A. We discuss the effect of non-stationary white
noise in App. B. All parameters and priors are listed in Tabs. 2 and 3.

###  4.2 Non-stationary QPO in white noise

We start by considering a combination of a QPO and white noise where the setup
is effectively identical to what we describe in 4.1, but we use
$S^{\mathrm{Q}}$ instead of $S^{\mathrm{RQ}}$ to omit red noise. We consider
both scenarios discussed in Sec. 3.1 where we either extend the time series
with zeros or with white noise. In terms of the parameter inference, we
deviate from the standard practice of using $S^{\mathrm{RWQ}}(f)$ as there is
no red noise component. Instead we compare $S^{\mathrm{WQ}}$ against the white
noise hypothesis $S^{\mathrm{W}}=C$.

First, we consider what happens when we extend the time series with white
noise. In principle, the results should follow our discussion in Sec. 3.1 for
QPOs in the white-noise-dominated part of the periodogram, i.e. the SNR should
decrease with $1/\sqrt{x}$ and the $\ln BF$ should decrease roughly inversely
to $x$. The results, which we display in the top panel of Fig. 2, are in
agreement with our expectations.

We now consider the case where we extend the time series with zeros, similar
to what we did to produce Fig. 1. As we discussed in Sec. 3, we expect the SNR
to increase with $\sqrt{x}$ and the $\ln BF$ roughly linearly. As we show in
the bottom panel of Fig. 2, this is almost perfectly the case for the $\ln BF$
and is also qualitatively true for the SNR.

Figure 2: SNR (blue) and $\ln BF$ (orange) vs the extension factor for a
simulated QPO plus white noise signal that is extended with white noise (top)
or zeros (bottom). The solid blue line displays the SNR of the maximum
likelihood point and the shaded blue region displays the SNR’s 90% credible
interval. In the top panel both SNR and $\ln BF$ decrease with increasing $x$
which corresponds to the QPO vanishing in noise. In the bottom panel the $\ln
BF$ increases almost perfectly linearly due to the increased number of
frequency bins that make up the QPO, whereas the SNR increases roughly with
the square root.

###  4.3 Non-stationary QPO in a simple transient

We investigate the effect of analyzing a non-stationary QPO that appears on
top of a deterministic transient flare shape in the presence of white noise.
The time series is displayed in the top panel of Fig. 3. The flare shape
mimics a red-noise continuum that overlaps with the QPO in the periodogram
(bottom panel of Fig. 3), though it is not well described by a Whittle
likelihood in the periodogram since it is not a stochastic process. We see
that there is conspicuously little fluctuation at low frequencies in Fig. 3.
For the flare shape we use a continuous exponential rise and fall described by

$\displaystyle
f(t)=A_{\mathrm{flare}}\begin{cases}\exp\left(\frac{t-t_{0}}{\tau_{\mathrm{r}}}\right)\,&\mathrm{if}\,t<t_{0}\\\
\exp\left(-\frac{t-t_{0}}{\tau_{\mathrm{f}}}\right)\,&\mathrm{if}\,t\geq
t_{0}\end{cases}$ (26)

where $A_{\mathrm{flare}}$ is the amplitude, $t_{0}$ is the peak time, and
$\tau_{\mathrm{r}}$ and $\tau_{\mathrm{f}}$ are the rise and fall timescales,
respectively. This model is similar to the FRED model which is popular to fit
time series of gamma-ray bursts (Norris et al., 2005). The advantage of using
an exponential shape for this simulation is that it produces a power law in
the periodogram which matches our red noise model. Since the flare is
deterministic, low frequencies in the periodogram are not well described by a
Whittle likelihood. To create this time series, we repeat the steps in Sec.
4.2, though we start with different parameters listed in Tab. 3, and add the
deterministic flare shape in the last step.

We perform the same test as in the previous section of extending the time
series starting from the $20\text{\,}\mathrm{s}$ segment in which the QPO is
present (inset in the top panel of Fig. 3). It is less intuitive why non-
stationarity bias may occur and how it will manifest in this scenario.
Extending beyond the peak of the flare might reduce the apparent significance
as it will provide the strongest contribution to the low-frequency noise
continuum. At the same time, we extend towards the tail of the flare, which
contributes little more than white noise. As we can see in Fig. 4, the SNR and
$\ln BF$ do indeed increase sharply until $x=3$, which corresponds to the
point where we extend the selected window to the peak of the transient. Past
this point low-frequency contributions from the rising edge of the transient
obscure the QPO in the periodogram and quickly lead to the QPO becoming
undetectable for $x>7$. We note that using a different set of parameters to
create the data can easily create situations where the $\ln BF$ continues to
increase past $x=3$.

Figure 3: Simulated data of a non-stationary QPO in a deterministic transient
flare shape in the presence of white noise. The time series (top) shows the
dominant transient and the QPO that is present for $20\text{\,}\mathrm{s}$ on
the tail of the flare (inset). We remove the deterministic trend of the flare
in the inset to make the QPO easily visible. The periodogram (bottom) shows a
visible QPO at $1\text{\,}\mathrm{Hz}$ and displays a low-frequency continuum
that arises due to the transient flare shape. The periodogram corresponds only
to the segment of the lightcurve that contains the QPO. Note that there is
much less variability than we would expect from the
$\chi^{2}_{2}$-distribution for the bins at low frequencies since the noise
continuum is not due to a stationary stochastic process but due to a
deterministic process.  Figure 4:  SNR (blue) and $\ln BF$ (orange) vs the
extension factor for the QPO shown in Fig. 3. The solid blue line displays the
SNR of the maximum likelihood point and the shaded blue region displays the
SNR’s 90% credible interval. Both the $\ln BF$ and the upper edge of the SNR’s
credible interval peak at $x=3$ (black dashed vertical line), which
corresponds to the point where the selected data segment exceeds past the peak
of the flare. The strong differences in $\ln BF$ and SNR across extension
factors shows that differences data selection can greatly influence the
results.

###  4.4 Non-stationary QPO in non-stationary red noise

In this scenario, we create a time series using a combination of red noise,
white noise, and a QPO using the $S^{\mathrm{RWQ}}$ model. We extend this time
series with either white noise or zeros. Fig. 5 displays the time series for
$|t|<$30\text{\,}\mathrm{s}$$ (top panel) and the periodogram as obtained for
$|t|<$10\text{\,}\mathrm{s}$$.

In the case of extending with zeros, we expect the same relation as in Sec.
4.2, which we recover in the bottom panel of Fig. 6. On the other hand, when
we extend with white noise the SNR first rises quickly until the extension
factor reaches around $x_{\mathrm{break}}=2$ (see Eq.24) and then turns
towards a slow descent proportional to $1/\sqrt{x}$. We note that the behavior
in the case of $x<x_{\mathrm{break}}$ for white noise is not well captured,
but greater extension factors would be less realistic and the simulated data
set aims to qualitatively show all effects.

We use models $S^{\mathrm{RW}}$ and $S^{\mathrm{RWQ}}$ as described earlier to
calculate posterior distributions, Bayes factors, and $\Delta BIC$s. As we
show in the bottom panel of Fig. 6, the $\ln BF$ increases almost perfectly
linearly when we extend with zeros. In the top panel, we see that if we extend
with white noise the $\ln BF$ also increases with $x$ even past
$x_{\mathrm{break}}$, though turns around eventually. The $\ln BF$ may be
increasing past $x_{\mathrm{break}}$ if the QPO feature in the periodogram is
better fit by the model for higher $x$ compared to the feature’s shape at
$x_{\mathrm{break}}$.

We have shown that overall we can see a non-stationarity bias if white noise
remains stationary, but red noise vanishes for part of the lightcurve. This is
in principle different from the effect in Sec. 4.3, as the red noise is due to
a stochastic rather than a deterministic process.

Figure 5: Simulated data of a non-stationary QPO and red noise in surrounding
stationary white noise. In the time series (top) we show how QPO and red noise
are transitioning smoothly into the surrounding part of the time series that
is pure white noise ($|t|>$10\text{\,}\mathrm{s}$$) due to the Hann window
that we have applied. In the periodogram (bottom) for
$|t|\leq$10\text{\,}\mathrm{s}$$ the QPO is clearly visible by eye.

Figure 6: SNR (blue) and $\ln BF$ (orange) vs the extension factor for
simulated data of a non-stationary QPO and red noise with white noise as
described in Sec. 4.4 for extension with white noise (top) and zeros (bottom).
The solid blue line displays the SNR of the maximum likelihood point and the
shaded blue region displays the SNR’s 90% credible interval.

###  4.5 Tests of stationarity

One noticeable property of the bias we have found is that it passes regular
tests employed to detect fitting validity. Nita et al. (2014) demonstrates a
$\chi^{2}$-like test for exponentially distributed data

$\chi^{2}_{\nu}=\frac{1}{\nu}\sum_{j=1}^{N/2}\left(1-\frac{I_{j}}{\hat{S}_{j}}\right)^{2}\,,$
(27)

where $\nu$ are the number of degrees freedom and $\hat{S}$ is the maximum
likelihood PSD. Analogously to the regular $\chi^{2}$-test,
$\chi^{2}_{\nu}\approx 1$ indicates an appropriate fit to the data, whereas
$\chi^{2}_{\nu}>1$ indicates underfitting, and $\chi^{2}_{\nu}<1$ overfitting.
Poor fits should in general occur with non-stationary behavior in the red
noise or QPO dominated part of the PSD. That is because neighboring bins will
be covariant and a fit may either align more closely or further away than
should be possible for independent bins, though $\chi^{2}\approx 1$ may still
randomly occur. However, in practice, non-stationarities will affect the low-
frequency part of the PSD while high frequencies are dominated by white noise.
Since frequency bins are linearly spaced, the vast majority of all bins are
almost always fitted well in terms of the $\chi^{2}$-statistic if the
periodogram reaches a white noise floor at high frequencies. As we show in
Fig. 7, the $\chi^{2}$-value for the entire PSD is very close to the
$\chi^{2}$-value for all the frequencies above the QPO frequency plus twice
its width. For comparison, setting $f_{0}=$1\text{\,}\mathrm{Hz}$$ with a
$40\text{\,}\mathrm{Hz}$ sampling frequency, as it is the case in the bottom
panel of Fig. 7, about $95\%$ off all frequency bins are greater than $f_{0}$.
Thus even if we overfit locally around the QPO, it is unlikely that this
greatly impacts the overall $\chi^{2}$-value.

We propose a local $\chi^{2}$-test, where we calculate the $\chi^{2}$-value
for just the frequencies within $2\sigma$ around the QPO frequency. We show
this test in Fig. 7 for data sets used in Secs. 4.2 (top) and 4.4 (bottom). We
note in the top panel that we obtain a reasonable $\chi^{2}\approx 0.92$ for
the local test if we just look at the stationary time series at $x=1$, and
obtain $\chi^{2}<0.6$ for large $x$. This method has the downside that it will
not work if there are only very few bins available as is the case for
$x\lesssim 6$ in the bottom panel of Fig. 7. We need at least seven data
points to compensate for the six degrees of freedom in the $S^{\mathrm{RWQ}}$
model. If more bins are available, this test can reliably flag whether the QPO
has been overfitted or underfitted. This test is required but not sufficient
to show that non-stationarities are not impacting the analysis. In general, we
may find $\chi^{2}\approx 1$ even for QPOs that are strongly affected by non-
stationarity bias. For example, a time series may in principle contain two
QPOs at similar frequencies and non-stationarity bias. In that instance, the
overfitting due to covariance of neighboring bins may be compensated by
underfitting of the underlying multiple QPO shapes and yield $\chi^{2}\approx
1$ for the local test.

Figure 7: Nita et al. (2014) $\chi^{2}$-tests in the case of a QPO and white
noise extended with zeros as described in Sec. 4.2 (top), and the case of QPO,
red noise and white noise extended with white noise as described in Sec. 4.4
(bottom). We see in the top panel that the overall fit (blue) and the fit for
frequencies above the QPO frequency (green) have $\chi^{2}\approx 1$,
indicating a good fit. The local $\chi^{2}$-test around the QPO frequency
$f_{0}$ indicates a good fit only for $x=1$ and falls much below that for
$x>1$, indicating overfitting. This implies that the local $\chi^{2}$-test can
detect non-stationarity bias in this instance. In the bottom panel, we display
some of the limitations of the local $\chi^{2}$-test when we deal with narrow
QPOs. For $x\lesssim 6$ the local $\chi^{2}$-test (orange) around $f_{0}$
swings widely as we have too few frequency bins within the range. Above the
$x\approx 6$ enough frequency bins accumulate within the range and we detect
the overfitting. We also note that the $\chi^{2}$-test for the entire PSD
(blue) and frequencies above the QPO frequency (green) indicate an overall
good fit despite us overfitting the QPO

##  5 Solar flare data

QPOs are a regularly reported phenomenon that have been observed in solar
flares over a wide range of wavelengths for decades, from radio waves to EUV
as well as soft and hard X-rays (see Nakariakov & Melnikov (2009); Van
Doorsselaere et al. (2016) for recent reviews). In this domain, they are
usually referred to as quasi-periodic pulsations (QPPs). Due to their
ubiquity, progress has been undertaken to automate their analysis using
Whittle likelihoods, e.g. with the Automated Flare Inference of Oscillations
(AFINO) method (Inglis et al., 2015, 2016; Hayes et al., 2020). An extensive
catalog of AFINO analyses is available
online333https://aringlis.github.io/AFINO/. AFINO defines models similar to
$S^{\mathrm{RW}}$, $S^{\mathrm{RWQ}}$, and $S^{\mathrm{BW}}$, although it uses
a Gaussian instead of a Lorentzian pulse as the QPO component, and a slightly
different broken power-law model. AFINO uses a computationally efficient Scipy
fitting routine to determine the maximum likelihood values and calculate the
$\Delta BIC$ with Eq. 8, though it does not calculate Bayes factors or
posterior distributions. The QPO model is “strongly favored” if $\Delta
BIC<-10$, relative both to the red noise and broken power-law model. This
ensures a small number of false-positive values, as has been demonstrated on
simulated data in Broomhall et al. (2019). Additionally, AFINO performs the
$\chi^{2}$-test for exponentially distributed data using Eq. 27. AFINO flags
the results if the implied p-value from the $\chi^{2}_{\nu}$ test is below or
above a set threshold.

For a comparative study on the impact of non-stationarity bias, it is
instructive to consider some of the most significant QPOs. This is because
these QPOs are visually identifiable from the lightcurve, which guides us in
their analysis. We look at a solar flare observed by the X-Ray Sensor (XRS)
onboard the GOES-15 satellite on May 12, 2013, at 20:17 UT for which AFINO
reports the highest significance for a QPO with $\Delta BIC=-451.6$ relative
to red noise and $\Delta BIC=-278.2$ relative to the broken power law. This
flare was of GOES magnitude M1.9 and originated from NOAA active region 11748,
which was located on the East limb of the Sun at the time. AFINO finds a high
quality fit for the QPO model with $\chi^{2}=1.06$ but flags the fits with the
broken power law and red noise due to their high $\chi^{2}$ values, likely due
to the very pronounced QPO present at $P=$12.6\text{\,}\mathrm{s}$$. We show
the 46 minutes long GOES x-ray lightcurve we analyze in Fig. 8. The QPO is
visible by eye on the tail end of the distribution (inset) and persists for
about $200\text{\,}\mathrm{s}$, much shorter than the
$2760\text{\,}\mathrm{s}$ data segment. This particular QPO may not be of
solar origin, since it is not observed in other flare-observing instruments.
However, a similar signature is coincidentally seen in the GOES magnetometer
data, suggesting a possible artifact. Regardless of the origin, the QPO is
clearly present and confined to a small portion of the data, providing an
ideal test case.

Figure 8: Solar flare x-ray lightcurve as observed with GOES in the
$1\text{\,}\mathrm{\SIUnitSymbolAngstrom}8\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
band. The figure shows the time selected as done with AFINO. The inset zooms
in on the tail of the lightcurve where the QPO is clearly visible.

To demonstrate the presence of an artificial amplification of the SNR, we
split the lightcurve into three segments. The first segment covers the time
before the QPO ($0\text{\,}\mathrm{s}1680\text{\,}\mathrm{s}$) (after 20:17
UT), the second selects just the QPO
($1680\text{\,}\mathrm{s}1880\text{\,}\mathrm{s}$), and the third segment
selects everything after the QPO
($1880\text{\,}\mathrm{s}2760\text{\,}\mathrm{s}$). We perform Bayesian
inference separately on each of these segments and the combined segment
independently with the model $S^{\mathrm{RW}}$ and $S^{\mathrm{RWQ}}$. We use
the same priors as in the studies on simulated data which we list in Tab. 2.
If the lightcurve were stationary, it would be valid to combine the Bayes
factors and $f_{0}$ posteriors as they represent independent draws from the
same distribution. However, as we show in Tab. 1, we are unable to detect the
reported QPO in the first and third segments. There is some weak evidence
towards a QPO in the first segment, though not at the reported
$12.6\text{\,}\mathrm{s}$ period. The third segment shows some weak support
for a QPO with $\ln BF=1.4$, but the maximum likelihood fit indicates that
this is rather due to a broad feature in the periodogram that is better fitted
with a broken power law. As we show in Tab. 1, both the first and third
segments are better fitted with $S^{\mathrm{BW}}$ than with
$S^{\mathrm{RWQ}}$, which indicates that it is likely no QPO present. In the
second segment we find the QPO independently with a very high significance
$\ln BF=27.5$, and it is clearly visible in the fitted periodogram in Fig. 9.
Finally, we analyze the combined segments together and find the QPO with $\ln
BF=229.4$, a significance much higher than in segment two. We note that we
find the QPO at a slightly different period of
$$12.34\text{\,}\mathrm{s}$\pm$0.06\text{\,}\mathrm{s}$$ compared to the value
reported of AFINO. This may be due to slightly different QPO modeling choices
between AFINO and our method.

We perform the local $\chi^{2}$-test we introduced in Sec. 4.5 for the entire
lightcurve and find $\chi^{2}=0.29$ for the frequency bins surrounding the
QPO. This indicates that the QPO has been overfitted and is non-stationary.

Segment |  $\ln BF$ $S^{\mathrm{RWQ}}$ |  $\ln BF$ $S^{\mathrm{BW}}$ |  $\Delta BIC$ $S^{\mathrm{RWQ}}$ |  $\Delta BIC$ $S^{\mathrm{BW}}$
---|---|---|---|---
$0\text{\,}\mathrm{s}1680\text{\,}\mathrm{s}$ | $3.0\pm 0.3$ | $6.8\pm 0.3$ | -11.3 | -5.0
$1680\text{\,}\mathrm{s}1880\text{\,}\mathrm{s}$ | $27.5\pm 0.2$ | $-0.7\pm 0.2$ | -66.8 | 7.8
$1880\text{\,}\mathrm{s}2760\text{\,}\mathrm{s}$ | $1.4\pm 0.3$ | $7.7\pm 0.2$ | -0.7 | -6.8
$0\text{\,}\mathrm{s}2760\text{\,}\mathrm{s}$ | $229.4\pm 0.3$ | $93.0\pm 0.3$ | -465.6 | -175.8
Table 1: Results from analyzing the selected segments of the solar flare
detected by GOES. All errors are given based on $1-\sigma$ confidence or
credible intervals. The segments are given in seconds after 20:17 UT. All $\ln
BF$ and $\Delta BIC$ values are calculated relative to $S^{\mathrm{RW}}$. We
reiterate that either a positive $\ln BF$ or a negative $\Delta BIC$ indicate
that $S^{\mathrm{RWQ}}$ or $S^{\mathrm{BW}}$ are preferred over
$S^{\mathrm{RW}}$. We find broadly that the $\Delta BIC$ values in the last
row are in agreement with what AFINO reported. Deviations are likely due to
differences, in our model and slightly different data selection. While in all
instances Bayes factors and $\Delta BIC$s give the same indication about the
preferred model, the significance differs to some extent. This may be in part
due to our wide prior choices.  Figure 9: Maximum likelihood fit of the
periodogram just for the second segment of the GOES lightcurve. The QPO is
clearly visible and highly significant. We draw the $2\sigma$ and $3\sigma$
detection limits based on the Bonferroni-corrected frequentist statistics. We
note that the QPO is located in the red noise dominated part of the PSD.

It is evident that the solar flare segment we have analyzed contains a QPO and
an overestimation of its significance is not as critical as it would be for a
marginal detection. Trying to find an instance where a marginal detection was
turned into a very confident detection due to non-stationarity bias would be
much harder. We would not be able to determine the location of the QPO in the
lightcurve by eye. Instead, we would have to take a systematic approach and
split the lightcurve into several smaller segments and determine if the QPO
exists in them individually. Of course, the significance in the individual
segments would always be lower, so it would be hard to determine whether this
is due to the reduced non-stationarity bias or because we split a persistent
QPO into multiple segments. As mentioned previously, we do not expect AFINO to
report many false detections, as has been established in Broomhall et al.
(2019). This is due to its conservative detection threshold $\Delta BIC<-10$.
Thus, if AFINO is coupled with a more careful data selection method it could
be possible to find more QPOs.

The non-stationarity bias in solar flare data arises most likely as a
combination of the effects we describe in Sec. 3, i.e. a non-stationary QPO,
non-stationary noise processes, and a deterministic overall flare shape all
contribute to some extent. We can not easily discern what the impact of each
of these effects is.

Non-automated analyses of QPOs in solar flares are also likely less affected
than the stretch of GOES data we analyze in this section. Typically, parts of
the lightcurve that are of interest are selected manually for detailed
analyses (e.g. Hayes et al., 2019), though this by no means provides a
guarantee that such results are without bias. As we have established
previously, the $\ln BF$ or the $\Delta BIC$ can grow linearly with $x$, thus
even a mild overselection of the lightcurve can lead to erroneous detections
and inferences.

##  6 Discussion and outlook

In this paper, we show that analyses based on Whittle likelihoods likely
overstates the significance of QPOs in lightcurves if the QPO or noise
features in time series are non-stationary. Specifically, this is the case if
data is included in the time series that contains only white noise when the
QPO is in the red-noise-dominated part of the PSD. The effect can also occur
if the QPO is in the white-noise-dominated part of the PSD, but the white
noise is non-stationary (e.g. if the noise arises due to photon counting) and
we include weak white noise from parts of the lightcurve that have low photon
counting rates. Selecting an appropriate time segment for the analysis of
transients is thus important but remains difficult to do generically and
without bias. A scientist’s natural intuition may be to select data
conservatively, i.e. select a longer time series while being aware that the
QPO may only be present for a part of it. The erroneous reasoning may be that
by selecting a longer segment they are not “cherry-picking” the segment with
the most pronounced QPO features. As we demonstrate in this paper, this choice
may increase the significance of the QPO artificially and lead to false
detections due to the non-stationarity bias.

There are some possible mitigation strategies that one may apply when the use
of periodograms is still preferred. However, this means that it is a lot
harder to devise a generic process with which to analyze time series. Firstly,
if we suspect that the QPO is indeed only present for part of the transient,
and the white noise level is relatively low compared to the low-frequency
continuum, as is indeed the case for x-ray lightcurves of solar flares, it is
reasonable to seek out the part of the lightcurve that looks by eye most
likely to contain a QPO, and then set the limits on the start and end of the
segment manually. However, this method requires significant human supervision
and thus will not scale to large data sets. It is open to other forms of
selection bias, too. Secondly, for such transients, we also suggest splitting
the lightcurve up into several parts to see if one can detect the QPO in all
segments or only in some of them. For the final analysis, one should only use
segments if there is some evidence for the QPO. Next, for shorter transients
where we suspect the QPO to be present for most of the time, we should
conservatively trim the lightcurve after the onset and before the end of the
segment. This way we are less likely to make the lightcurve non-stationary.
Finally, we have found that regular $\chi^{2}$-like tests are unsuitable to
detect non-stationary bias since only a few bins around the QPO are affected.
The overall $\chi^{2}$-value is mostly determined by the far greater share of
frequency bins that lie in the white-noise-dominated part of the PSD. We
outline that alternatively calculating $\chi^{2}$ based solely on the
frequencies surrounding the QPO can detect overfitting, which is a hallmark of
statistically non-independent frequency bins.

Aside from periodograms, other methods that can be used for the search for
QPOs. Wavelet transforms are a popular way to analyze the time series by
convolving a wavelet function with a time series of interest (Torrence &
Compo, 1998). Wavelet analyses are not restricted to stationary data sets and
may help us to find the specific times when a QPO occurs. However, the
statistics of detecting QPOs with wavelets in the presence of red noise
remains contentious since detrending methods are likely to lead to false
detections (Auchère et al., 2016). Alternatively, we can avoid frequency-
domain methods altogether and model the time series as a Gaussian process with
some mean function. The main drawback of Gaussian processes is that they have
computational complexity $\mathcal{O}(N^{3})$ and thus are only suitable for
short time series in the general case (Rasmussen & Williams, 2006), though
progress has been made to reduce complexity to $\mathcal{O}(N)$ for
stationary, complex exponential kernels or combinations thereof (Foreman-
Mackey et al., 2017). We note that we expect the same non-stationarity bias if
we apply a stationary Gaussian process kernel to a non-stationary time series
since the Whittle likelihood explicitly derives from a Gaussian process
likelihood. We have also found this result empirically in some preliminary
analyses (Huebner et al. 2021, in prep).

One advantage of Gaussian processes is that they allow us to fit the overall
shape of the transient and the stochastic process simultaneously, instead of
lumping the shape of the transient into the red noise, which can help us to
prevent the bias demonstrated in Sec. 4.3. It thereby allows us to avoid
nonparametric methods that are likely to create false periodicities (Auchère
et al., 2016). Additionally, the Gaussian process likelihood can explicitly
take in the known variance for flux values within a lightcurve if they are
available. Thus, they also provide a natural resolution to the bias shown in
App. B. Finally, Gaussian processes allow us to model relatively simple, non-
stationary extensions to the fast, stationary models that are popular now
(Huebner et al. 2021, in prep). Further development of Gaussian process
methods may eventually allow us to create more complex models of QPOs and
their possible intermittency within transients.

#### Acknowledgments

This work is supported through ARC Centre of Excellence CE170100004. P.D.L. is
supported through the Australian Research Council (ARC) Future Fellowships
FT160100112, ARC Discovery Project DP180103155. D.H. is supported by the Women
In Science Excel (WISE) programme of the Netherlands Organisation for
Scientific Research (NWO). This work was performed on the OzSTAR national
facility at Swinburne University of Technology. The OzSTAR program receives
funding in part from the Astronomy National Collaborative Research
Infrastructure Strategy (NCRIS) allocation provided by the Australian
Government.

##  Appendix A Priors and simulated data parameter tables

We list the priors and parameters used for all the studies in Secs. 4 and 5 in
Tabs. 2 and 3.

Parameter | Description | Prior | Minimum | Maximum
---|---|---|---|---
$A$ | Red noise amplitude | LogUniform | $\exp(-30)$ | $\exp(30)$
$\alpha$, $\alpha_{1}$, $\alpha_{2}$ | red noise / BPL spectral index | Uniform | $0$ | $10$
$B$ | QPO amplitude | LogUniform | $\exp(-60)$ | $\exp(60)$
$f_{0}$ | QPO central freq. | LogUniform | $2\Delta f$ | $f_{\max}$
$\sigma$ | QPO HWHM | ConditionalLogUniform | $\Delta f/\pi$ | $0.25f_{0}$
$C$ | White noise amplitude | LogUniform | $\exp(-30)$ | $\exp(30)$
$\delta$ | Breaking freq. for BPL | LogUniform | $\Delta f$ | $f_{\max}$
Table 2: Parameters and priors used throughout all studies in this paper. Priors are referred to by their implementation in Bilby. We select wide priors in $A$, $B$, and $C$ for simplicity. Although the selection of priors this wide may have a slight impact on the Bayes factor calculation, it will not qualitatively change the scaling of Bayes factor with $x$. We also set a wide prior on $\alpha$, or $\alpha_{1}$ and $\alpha_{2}$, and enforce $\alpha_{2}<\alpha_{1}$ using a Constraint prior in Bilby to avoid degeneracies. The prior on $f_{0}$ is motivated by the available frequencies in the periodogram. We set the minimum $f_{0}$ at twice the difference between neighboring frequencies $\Delta f$ to ensure better convergence, and refer to the highest frequency in the periodogram as $f_{\max}$. We set $\sigma$, the width of the QPO at half width half maximum (HWHM), to be conditional uniform in log between $\Delta f/pi$, i.e. on the scale of a single frequency bin, and $0.25f_{0}$, to prevent it from converging towards wider features in the periodogram. For the broken power law analysis of solar flare data, we use priors for $\alpha_{1,2}$ identical to $\alpha$, and the listed $\delta$ prior Parameter | Sec. 4.2 zeros | Sec. 4.2 white noise | Sec. 4.3 | Sec. 4.4 | App. B
---|---|---|---|---|---
$A$ | - | - | - | 4 | -
$\alpha$ | - | - | - | 2 | -
$B$ | 3 | 15 | 80 | 8 | 100
$f_{0}$ | 5 | 5 | 1 | 1 | 5
$\sigma$ | $0.1\text{\,}\mathrm{Hz}$ | $0.1\text{\,}\mathrm{Hz}$ | $0.01\text{\,}\mathrm{Hz}$ | $0.02\text{\,}\mathrm{Hz}$ | $0.1\text{\,}\mathrm{Hz}$
$C$ | 2 | 2 | 2 | 2 | -
$A_{\mathrm{flare}}$ | - | - | 20000 | - | $1000\text{\,}{\mathrm{s}}^{-1}$
$t_{0}$ | - | - | $70\text{\,}\mathrm{s}$ | - | $200\text{\,}\mathrm{s}$
$\tau_{\mathrm{r}}$ | - | - | $10\text{\,}\mathrm{s}$ | - | -
$\tau_{\mathrm{f}}$ | - | - | $20\text{\,}\mathrm{s}$ | - | -
$\sigma_{\mathrm{flare}}$ | - | - | - | - | $20\text{\,}\mathrm{s}$
$c_{0}$ | - | - | - | - | $10\text{\,}{\mathrm{s}}^{-1}$
Segment length | $400\text{\,}\mathrm{s}$ | $400\text{\,}\mathrm{s}$ | $200\text{\,}\mathrm{s}$ | $400\text{\,}\mathrm{s}$ | $400\text{\,}\mathrm{s}$
Sampling frequency | $40\text{\,}\mathrm{Hz}$ | $40\text{\,}\mathrm{Hz}$ | $40\text{\,}\mathrm{Hz}$ | $40\text{\,}\mathrm{Hz}$ | $40\text{\,}\mathrm{Hz}$
Table 3: Values used to create the simulated data in Sec. 4. The column heads
refer to the subsection in which this set of parameters was used. We use two
different sets of parameters in Sec. 4.2 depending on whether we extend using
zeros or more white noise.

##  Appendix B Deterministic processes impact white noise

White noise observed in astrophysical lightcurves using photon counting does
not arise due to intrinsic properties of the source, but rather due to the
Poissonian noise nature of photon counting. Concretely, given a rate
$\lambda$, the distribution of the number of observed photons $k$ in a unit of
time is Poissonian

$\mathrm{Pois}(k;\lambda)=\frac{\lambda^{k}e^{-\lambda}}{k!}\,.$ (B1)

Photon counting noise thus scales proportionally to the standard deviation of
the Poisson distribution $\sqrt{\lambda}$. This relation implies that photon-
counting noise throughout a transient does not remain constant. We demonstrate
a limit of this in Sec. 3.1 when we consider what happens if we extend a time
series with zeros. This limit corresponds to the case in which a detector sees
no photons at all. In realistic detectors, we will not reach this limit as
there are at least some background photons, though we can get close enough to
it for it to lead to wrong inferences about QPO significance.

The above effect points to a deficiency in periodograms more generally. For
the calculation of a periodogram we use the values $x(t_{i})$ of a discrete
time series, but we do not include their variances $\Delta x(t_{i})$ if they
are available. In effect, this implies a loss of information as the white
noise is already encoded in the variance associated with the photon counts.
Instead, we infer it after the fact via power spectral density estimation
independently and lose knowledge about non-stationarity in the white noise.

To demonstrate this effect, we construct a transient similar to Sec. 4.3,
though in a way that its low-frequency contributions do not overlap with the
QPO frequency, as to demonstrate that this effect arises due to non-
stationarity in the white noise. Concretely, we construct the transient by
generating a Gaussian profile plus a constant

$f(t;A_{\mathrm{flare}},t_{0},\sigma_{\mathrm{flare}})=A_{\mathrm{flare}}\exp\left(-\frac{(t-t_{0})^{2}}{2\sigma_{\mathrm{flare}}^{2}}\right)+c_{0}\,,$
(B2)

where $A_{\mathrm{flare}}$ is the amplitude, $t_{0}$ is the time of the peak,
$\sigma_{\mathrm{flare}}$ is the width, and the constant $c_{0}$ represents a
possible background count rate. We add a $20\text{\,}\mathrm{s}$ non-
stationary QPO in the same manner as we describe in Sec. 4.3. We create photon
counts, which we show in Fig. 10, by simulating the Poisson process using the
implementation in scipy.stats.poisson.rvs. The specific parameters of this
time series are listed in Tab. 3. By construction, the QPO is set on the top
of the Gaussian profile, which corresponds to the highest level of photon
counting noise. The Gaussian profile, unlike the exponential from Sec. 4.3,
contributes powers at lower frequencies than the exponential, roughly up to
$1/\sigma_{\mathrm{flare}}\approx$0.05\text{\,}\mathrm{Hz}$$. Additionally, we
set our QPO frequency at $5\text{\,}\mathrm{Hz}$ and cut off frequencies below
$0.5\text{\,}\mathrm{Hz}$ so that we can avoid using a red noise component in
our modeling in the same way we did in Sec. 4.2.

On the tails of the flare, we expect the same effect as we observe when we
extend with zeros in our other simulations, though Fig. 11 clearly shows that
we artificially increase the significance already with very small extension
factors. This is possible since the flare quickly falls off on either side and
the added photon counting noise does not compensate for the increased number
of frequency bins that contain the QPO. Overall, this effect would be hard to
account for just based on better methods for setting data cuts because varying
levels of white noise are expected in many transient lightcurves.

Figure 10: Simulated data of a non-stationary QPO in a deterministic transient
flare shape in the presence of white photon counting noise. The time series
(top) shows the dominant transient and the QPO that is present for
$20\text{\,}\mathrm{s}$ on the top of the flare (inset). The periodogram for
the entire $400\text{\,}\mathrm{s}$ time series (bottom) has a relatively wide
($\sigma=0.1$) QPO at $5\text{\,}\mathrm{Hz}$ (barely visible). For our
analysis, we cut off all frequencies below $0.5\text{\,}\mathrm{Hz}$ to avoid
any effects arising from the power in low frequencies due to the flare shape.
Figure 11: We calculate the Bayes factors for the presence of a QPO in a
transient flare for different extension factors with $x=1$ corresponding the
$20\text{\,}\mathrm{s}$ inset from the top of Fig. 10. While we already start
with very strong evidence of a QPO at $x=1$, the Bayes factors continue to
increase as we extend the time series to include more of the transient.

## References

  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, The Astrophysical Journal Supplement Series, 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Auchère et al. (2016) Auchère, F., Froment, C., Bocchialini, K., Buchlin, E., & Solomon, J. 2016, The Astrophysical Journal, 825, 110, doi: 10.3847/0004-637X/825/2/110
  * Barret & Vaughan (2012) Barret, D., & Vaughan, S. 2012, 131, doi: 10.1088/0004-637X/746/2/131
  * Bonferroni (1936) Bonferroni, C. 1936, Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze, 8, 3
  * Broomhall et al. (2019) Broomhall, A. M., Davenport, J. R., Hayes, L. A., et al. 2019, arXiv, doi: 10.3847/1538-4365/ab40b3
  * Cenko et al. (2010) Cenko, S. B., Butler, N. R., Ofek, E. O., et al. 2010, Astronomical Journal, 140, 224, doi: 10.1088/0004-6256/140/1/224
  * Chatfield (1975) Chatfield, C. 1975, The Analysis of Time Series: Theory and Practice. (Chapman & Hall ; distributed by Halsted Press)
  * Covino et al. (2020) Covino, S., Landoni, M., Sandrinelli, A., & Treves, A. 2020, The Astrophysical Journal, 895, 122, doi: 10.3847/1538-4357/ab8bd4
  * De Luca et al. (2010) De Luca, A., Esposito, P., Israel, G. L., et al. 2010, Monthly Notices of the Royal Astronomical Society, 402, 1870, doi: 10.1111/j.1365-2966.2009.16012.x
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, The Astronomical Journal, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Foster (1996) Foster, G. 1996, The Astronomical Journal, 112, 1709, doi: 10.1017/CBO9781107415324.004
  * Gelman et al. (2013) Gelman, A., Carlin, J. B., Stern, H. S., et al. 2013, 677.  http://files/25/Gelmanm.fl.-BayesianDataAnalysisThirdedition(witherrors.pdf
  * Hayes et al. (2019) Hayes, L. A., Gallagher, P. T., Dennis, B. R., et al. 2019, arXiv, doi: 10.3847/1538-4357/ab0ca3
  * Hayes et al. (2020) Hayes, L. A., Inglis, A. R., Christe, S., Dennis, B., & Gallagher, P. T. 2020, The Astrophysical Journal, 895, 50, doi: 10.3847/1538-4357/ab8d40
  * Huppenkothen et al. (2014a) Huppenkothen, D., Watts, A. L., & Levin, Y. 2014a, Astrophysical Journal, 793, doi: 10.1088/0004-637X/793/2/129
  * Huppenkothen et al. (2012) Huppenkothen, D., Watts, A. L., Uttley, P., et al. 2012, Astrophysical Journal, 768, doi: 10.1088/0004-637X/768/1/87
  * Huppenkothen et al. (2014b) Huppenkothen, D., D’Angelo, C., Watts, A. L., et al. 2014b, Astrophysical Journal, 787, doi: 10.1088/0004-637X/787/2/128
  * Huppenkothen et al. (2017) Huppenkothen, D., Younes, G., Ingram, A., et al. 2017, The Astrophysical Journal, 834, 90, doi: 10.3847/1538-4357/834/1/90
  * Inglis et al. (2016) Inglis, A. R., Ireland, J., Dennis, B. R., Hayes, L., & Gallagher, P. 2016, The Astrophysical Journal, 833, 284, doi: 10.3847/1538-4357/833/2/284
  * Inglis et al. (2015) Inglis, A. R., Ireland, J., & Dominique, M. 2015, The Astrophysical Journal, 798, 108, doi: 10.1088/0004-637X/798/2/108
  * Ingram & Motta (2019) Ingram, A. R., & Motta, S. E. 2019, New Astronomy Reviews, 85, 101524, doi: 10.1016/j.newar.2020.101524
  * Ireland et al. (2015) Ireland, J., McAteer, R. T., & Inglis, A. R. 2015, Astrophysical Journal, 798, doi: 10.1088/0004-637X/798/1/1
  * Israel et al. (2005) Israel, G. L., Belloni, T., Stella, L., et al. 2005, The Astrophysical Journal, 628, L53, doi: 10.1086/432615
  * Iwakiri et al. (2010) Iwakiri, W., Ohno, M., Kamae, T., et al. 2010, AIP Conference Proceedings, 1279, 89, doi: 10.1063/1.3509358
  * Lachowicz & Done (2010) Lachowicz, P., & Done, C. 2010, Astronomy and Astrophysics, 515, A65, doi: 10.1051/0004-6361/200913144
  * Miller et al. (2019) Miller, M. C., Chirenti, C., & Strohmayer, T. E. 2019, The Astrophysical Journal, 871, 95, doi: 10.3847/1538-4357/aaf5ce
  * Morris et al. (2010) Morris, D., Battista, F., Dhuga, K., & MacLachlan, G. 2010, AIP Conference Proceedings, 1279, 394, doi: 10.1063/1.3509323
  * Nakariakov & Melnikov (2009) Nakariakov, V. M., & Melnikov, V. F. 2009, Space Science Reviews, 149, 119, doi: 10.1007/s11214-009-9536-3
  * Nita et al. (2014) Nita, G. M., Fleishman, G. D., Gary, D. E., Marin, W., & Boone, K. 2014, Astrophysical Journal, 789, doi: 10.1088/0004-637X/789/2/152
  * Norris et al. (2005) Norris, J. P., Bonnell, J. T., Kazanas, D., et al. 2005, The Astrophysical Journal, 627, 324, doi: 10.1086/430294
  * Rasmussen & Williams (2006) Rasmussen, C., & Williams, C. 2006, Gaussian Processes for Machine Learning, Vol. 7 (MIT Press)
  * Romero-Shaw et al. (2020) Romero-Shaw, I. M., Talbot, C., Biscoveanu, S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 3295, doi: 10.1093/mnras/staa2850
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833, doi: 10.1214/06-BA127
  * Speagle (2019) Speagle, J. S. 2019, arXiv, 3158, 3132, doi: 10.1093/mnras/staa278
  * Strohmayer & Watts (2005) Strohmayer, T. E., & Watts, A. L. 2005, The Astrophysical Journal, 632, L111, doi: 10.1086/497911
  * Strohmayer & Watts (2006a) —. 2006a, The Astrophysical Journal, 653, 593, doi: 10.1086/508703
  * Strohmayer & Watts (2006b) —. 2006b, The Astrophysical Journal, 653, 593, doi: 10.1086/508703
  * Tarnopolski & Marchenko (2021) Tarnopolski, M., & Marchenko, V. 2021, The Astrophysical Journal, 911, 20, doi: 10.3847/1538-4357/abe5b1
  * Timmer & König (1995) Timmer, J., & König, M. 1995, Astronomy and Astrophysics, 300, 707
  * Torrence & Compo (1998) Torrence, C., & Compo, G. P. 1998, Bulletin of the American Meteorological Society, 79, 61, doi: 10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2
  * van der Klis (1989) van der Klis, M. 1989, Annual Review of Astronomy and Astrophysics, 27, 517, doi: 10.1146/annurev.aa.27.090189.002505
  * Van Doorsselaere et al. (2016) Van Doorsselaere, T., Kupriyanova, E. G., & Yuan, D. 2016, Solar Physics, 291, 3143, doi: 10.1007/s11207-016-0977-z
  * Vaughan (2005) Vaughan, S. 2005, Astronomy and Astrophysics, 431, 391, doi: 10.1051/0004-6361:20041453
  * Vaughan (2009) —. 2009, Monthly Notices of the Royal Astronomical Society, 402, 307, doi: 10.1111/j.1365-2966.2009.15868.x
  * Watts & Strohmayer (2006) Watts, A. L., & Strohmayer, T. E. 2006, The Astrophysical Journal, 637, L117, doi: 10.1086/500735
  * Weakliem (1999) Weakliem, D. L. 1999, Sociological Methods and Research, 27, 359, doi: 10.1177/0049124199027003002
  * Whittle (1951) Whittle, P. 1951, Hypothesis testing in time series analysis (Almqvist & Wiksell)
  * Zhu & Thrane (2020) Zhu, X. J., & Thrane, E. 2020, arXiv, doi: 10.3847/1538-4357/abac5a
  * Zimovets et al. (2021) Zimovets, I. V., McLaughlin, J. A., Srivastava, A. K., et al. 2021, Space Science Reviews, 217, 66, doi: 10.1007/s11214-021-00840-9


# Open-Source Memory Compiler for Automatic RRAM Generation and Verification

Dimitrios (Dimitris) Antoniadis∗, Peilong Feng∗†, Andrea Mifsud∗†, Timothy G.
Constandinou∗†‡ ∗Department of Electrical and Electronic Engineering, Imperial
College London, SW7 2BT, UK
†Centre for Bio-Inspired Technology, Institute of Biomedical Engineering,
Imperial College London, SW7 2AZ, UK
‡Care Research & Technology Centre, UK Dementia Research Institute, UK
{dimitris.antoniadis20, peilong.feng14, a.mifsud,
<EMAIL_ADDRESS>

###### Abstract

The lack of open-source memory compilers in academia typically causes
significant delays in research and design implementations. This paper presents
an open-source memory compiler that is directly integrated within the Cadence
Virtuoso environment using physical verification tools provided by Mentor
Graphics (Calibre). It facilitates the entire memory generation process from
netlist generation to layout implementation, and physical implementation
verification. To the best of our knowledge, this is the first open-source
memory compiler that has been developed specifically to automate Resistive
Random Access Memory (RRAM) generation. RRAM holds the promise of achieving
high speed, high density and non-volatility. A novel RRAM architecture,
additionally is proposed, and a number of generated RRAM arrays are evaluated
to identify their worst case control line parasitics and worst case settling
time across the memristors of their cells. The total capacitance of lines SEL,
N and P is 5.83 fF/cell, 3.31 fF/cell and 2.48 fF/cell respectively, while the
total calculated resistance for SEL is 1.28 /cell and 0.14 /cell for both N
and P lines.

##  I Introduction

Random Access Memories (RAMs) and processing units are critical components of
a digital system. Commonly known volatile memories are Static Random Access
Memories (SRAMs) and Dynamic Random Access Memories (DRAMs) [1, 2]. SRAMs are
of particular interest since they can be used in System-On-Chip (SoC),
Application Specific Integrated Circuit (ASIC) and microprocessor designs [3,
4].

Memory design is however a tedious task for an Integrated Circuit (IC) design
engineer, and requires plethora of time and multiple design cycles [3, 4].
Additionally, commercial Process Design Kits (PDKs) that are available to
academia or small and medium size enterprises (SMEs), do not provide memory
compilers. Access to memory IP is typically licensed through third party
vendors, often providing just limited flexibility and reconfigurability [3].
Modern applications require customised memories (such as single port, dual
port, multi port) [5]. Depending on the desired performance, different
fabrication technology may be required [3], therefore, memory compilers should
be able to create scalable and customisable memories [3, 4]. Driven by the
above mentioned needs, various attempts have been made to create memory
compilers [3, 4, 5, 6].

Conventionally, memory compilers mainly deal with SRAMs or other volatile
memories [3, 4, 5, 6, 7, 8]. These memories can have advantages such as high
packing density and high read/write speed. Nevertheless, the volatility of the
states of their cells produces the necessity of constantly refreshing their
cell values (DRAMs) or restoring them on power-up (SRAMs). As a result the
energy consumption is high [1, 9]. On the other hand, commercially available
non-volatile memories (flash memories), suffer of low speed, high write
voltage and low endurance [1]. Both volatile and non-volatile memories face
the problem of the memory wall, as the write/read speed of memories does not
scale at the same rate as the operation frequency of processing units does
[10, 11, 2].

It is clear that novel material memories exploiting new memory architectures
are necessary to achieve non-volatility, scalability, high speed, high density
and low power performance [1, 2, 11, 12, 13]. These devices can be
characterised as emerging Non-Volatile Memory (NVM) technologies [1, 12]. Such
a device is the memristor [13, 14, 9]. A promising type of memory using
memristors and addressing the above mentioned problems is the Resistive RAM
(RRAM) [2, 9, 15].

Figure 1: Architecture of proposed RRAM: Top Right - RRAM array. The size of
the array is $M=2^{Y}$ rows and $N=b2^{X}$ columns. The blue coloured cells
represent a word of $b$ bits. Left - A $Y$ bits to $2^{Y}$ bits Decoder to
choose the address (enable horizontal SEL line) of the desired word. Bottom -
two $b2^{X}$ to $b$ De/Multiplexers (enable vertical P, N lines) to provide
the necessary connections based on the read/write operation. These are
connected to the Sense Amplifier Array and the Write Circuits Array shown in
orange and purple colour respectively.

Even though, a RRAM compiler has been presented before [15], the authors of
this paper acknowledge the need of an open-source RRAM compiler for academia
and industry. The compiler should be able to generate a RRAM and its
peripheral circuits. Additionally, the compiler should provide automatic
layout verification and memory characterisation. The proposed tool should be
independent of fabrication technology and should provide reconfigurability
options on the architectures and the circuits of the memory. Such a tool is
essential in order to investigate novel properties and capabilities of RRAMs.
This is developed as part of the Functional Oxide Reconfigurable Technologies
(FORTE) programme 111For latest news on FORTE and resources visit
www.forte.ac.uk. FORTE’s aims are threefold: develop and optimise the
manufacturing processes of RRAM devices; create the technology tools and
design rules required for integration with CMOS; and showcase potential
applications of CMOS/Memristor integrated circuits [9].

This paper presents an early version of the first open-source RRAM compiler
which automatically creates and verifies the schematic and the layout of a M
$\times$ N dimensions RRAM. Additionally, it extracts parasitics and generates
a post layout view. Suitable files with parasitic information are generated
which can be imported in MATLAB and further information on control lines
capacitance and resistance can be obtained. The RRAM compiler has been
published on Github:

https://github.com/akdimitri/RRAM_COMPILER

The proposed architecture of the RRAM of the compiler is shown on Fig. 1. The
subcircuits are quite similar to an SRAM [16], however, RRAM architecture
requires additional peripheral circuitry to control the RRAM cells which
operate on very strict voltage conditions as explained in Section II. Section
III describes the operation flow of the proposed compiler. Section IV compares
the characteristics of the various RRAM macros created by the compiler.
Finally, Section V discusses future improvements, challenges and concludes the
results.

##  II RRAM Architecture Overview

The proposed RRAM cell consists of a memristor and an NMOS transistor as shown
in Fig. 3. In order to write a value to the cell, the SEL line has to be set
to high voltage and an appropriate voltage difference has to be applied
between P and N terminals (depending on memristor characteristics). In order
to read the value of the memristor, SEL line has to be set to high voltage and
an extremely small voltage difference has to be applied between P and N
terminals, ensuring the voltage difference is not large enough to change the
resistive state of the memristor. Thus, it is made clear that suitable
peripheral circuits have to be designed in order to accurately control these
sensitive pins.

Figure 2: Read and Write operations. Assuming an unknown Resistive State RS on
the memristor, at $t=1$ a write signal appears on write/read (W/R) port and
Low Resistance ($L$) is received on Input/Output (I/O) port. At $t=3$ the read
operation of the cell begins. Note the small voltage applied on P line while
reading function of RRAM. At $t=8$, a high voltage pulse is applied on SEL and
N line in order to write a High Resistance ($H$) value.

The RRAM cells can form a RRAM array by sharing horizontal SEL lines and
vertical P and N lines as it is shown on Fig. 3. Based on this topology a M
$\times$ N RRAM array can be generated. Given that a word length is $b$ bits,
a $2^{Y}\times b2^{X}$ RRAM array can be generated, where $M=2^{Y}$ and
$N=b2^{X}$ with resemblance to the arrays of Figs. 1 and 3.

Figure 3: RRAM array $4\times 4$. RRAM rows share SEL line while RRAM columns
share P (blue) and N (red) lines. One memory cell is shown on the lower left
part of the image surrounded by dashed line.

Thus for both read and write operation, a SEL line decoder is needed to select
the desired horizontal line and multiplexers to connect the P and N vertical
lines of the desired word memory cells to read/write circuits. Since the array
has $2^{Y}$ SEL lines, a word of $Y$ bits is enough to describe the address of
the desired SEL line. Therefore a decoder of $Y$ to $2^{Y}$ is required.

Regarding P and N lines, it is essential to take under consideration the read
and write operations. During the read operation, a small voltage (less than
0.5 V) has to be applied on the desired P lines. Therefore, a multiplexer
$b2^{X}$ to $b$ is necessary to connect the selected memory word to read
circuits. A word of $X$ bits is enough to control the multiplexer. Since the
voltage difference has to be extremely small, the N pins can be connected to
ground (See Fig. 1 ground on N Multiplexer). Therefore, a signal can be used
(not shown on figure) during read operation to ground all N pins. On the other
hand, during write operation a larger voltage value has to be applied across P
and N pins, thus in this case both desired P and N lines will be connected to
write circuits using $b2^{X}$ to $b$ multiplexers. A simplified diagram of
Read and Write Operations is shown on Fig. 2. In Fig. 1 N multiplexer is
connected with pink line with the write circuits as it is needed only for
write operation. On the other hand, P multiplexer is connected on both read
and write circuits, thus a black line is used on the figure. The write
circuits are shown on Fig. 1 with WR blocks.

As far as the read operation is concerned, a novel sense amplifier has to be
designed. In contrast to conventional SRAM sense amplifiers and other sense
amplifiers for NVMs, which apply a minimum voltage of $1/3$_VDD_ or greater on
the bit line [16, 17, 18], the sense amplifier of the proposed RRAM has to
maintain an extremely small voltage on P line throughout the whole read
operation, in order not to change the resistive state of the memristor.
Potentially, no precharge circuit is required in contrast to SRAM
implementations [3]. The sense amplifiers array is shown on Fig. 1 with SA
blocks. These operations require a control logic circuit that will supervise
the whole operation of the RRAM circuits. Besides the abstract system
architecture presented on this section, further circuits are required, such as
Tri-state input/output ports. These kind of circuits and additional peripheral
circuits are essential for the RRAM operation and will be investigated and
added at a later stage of the Compiler design.

Figure 4: Generated $128\times 128$ RRAM array layout. On the top right of the
picture, a zoomed area is shown (poly, metal 1, metal 4). It depicts a memory
cell (1T1R) and the shared connections.

##  III Compiler Analysis

The proposed compiler was created using SKILL language which can be used to
customise and extend Cadence Design Environment (IC6.1.8)[19, 20]. In this
case, multiple SKILL scripts and functions have been created and integrated in
a main function SKILL script which performs the automatic memory generation
and layout verification. Both the main function and subfunctions are generic
functions, highly customisable and reconfigurable.

The arguments of the compiler are the $M$ x $N$ dimensions and the name of the
library where the generated RRAM should be saved. Additionally, the database
representation [19] of the RRAM cell and the Calibre cell-map file can be
provided as arguments or they can be directly hard-coded.

A simplified compilation flow is shown on Fig. 5. Initially the compiler
creates the schematic, the layout and the symbol of a row of RRAM Cells of
size $1\times N$. Based on the row cellview, the RRAM macro is generated.

Subsequently, the compiler creates three folders inside the cellview folder
for DRC, LVS and PEX respectively. The necessary files (runset files, Spice
netlist, GDS) are generated by the compiler in order to invoke Calibre DRC,
LVS and PEX in batch mode through the Command Interpreter Window (CIW) of
Virtuoso [21, 22]. After these operations are finished, the compiler invokes
Calibre View Setup in batch mode and creates a Calibre post-layout view of the
RRAM. The parasitics output files of the compiler are located inside the
cellview folders.

Figure 5: Flowchart of RRAM compiler operations.

The layout (metals and poly layers) of a generated $128\times 128$ RRAM array
is presented on Fig. 4. The layout size is 294.42 µm$\times$642.41 µm.

##  IV Results

Using the compiler, a number of $N$ x $N$ square RRAM arrays were created. By
importing the output files of the compiler in MATLAB, the worst case
Capacitances and Cross Capacitances ($C+CC$) were extracted. The results are
shown on Fig. 6a. Fig. 6a shows that $C+CC$ rises linearly with respect to the
number of cells of a line. The SEL lines have greater capacitance. The memory
cell has greater horizontal length compared to its vertical width with
resemblance to the memory cell layout of Fig. 4. Lines N and P are shared
vertically and occupy less area, thus they suffer of less capacitance. It is
worth mentioning that the SEL line capacitance does not include the additional
capacitance of the gates of the transistor switches and also that the
additional capacitance due to the connections of the cell to the memristor has
not been taken into account.

Figure 6: (a) Worst case Capacitance and Cross Capacitance ($C+CC$) of lines
SEL, P, N with respect to the number of Row/Column cells of RRAM. The increase
rates are 5.83 fF/cell, 3.31 fF/cell and 2.48 fF/cell for SEL, N and P lines
respectively. (b) Total calculated resistance of SEL, N, P lines with respect
to the number of cells of the line. The increase rate is 1.28 /cell for SEL
line and 0.14 /cell for N and P lines. (c) Worst case settling time of voltage
across the memristor of a cell. This case is for $V_{N}=$ _VDD_ , $V_{P}=$
_GND_ and SS model library. The curve fits the exponential curve $5.223\times
10^{-10}exp(0.004207x)$.

Additionally, the worst case scenario settling time of voltage across the
memristor on the post layout views of RRAM cells was investigated. This case
was found for _VDD_ on N pin and _VSS_ on P pin. As it is expected, the worst
case was found for the maximum high resistive state where R=5 M. Based on
these conditions a testbench was set in Virtuoso and a corner analysis was
conducted for a number of RRAM arrays. The settling time range was set equal
to 1$\%$ of the final value. Slowest settling time was achieved for Slow NMOS,
Slow PMOS (SS) model library. Results are shown on Fig. 6c. For small array
row/column size the settling time is almost stable around 550 ps, however as
the size grows the settling time rises exponentially. This is because both the
resistance and the capacitance of the lines increase and therefore the loading
time of the line increases.

Based on the SPICE model of the technology library for layout layers, the
resistance of lines SEL, N and P were calculated. N and P lines have the same
layout shape, thus it was found that they have the same resistance. The
results are shown on Fig. 6b. The resistance also rises linearly with respect
to the number of cells of the line.

The RRAM cell used for the results produces arrays of density equal to 0.082
Mb/mm2 for a 180 nm technology. The RRAM cells have approximately the same
size with SRAM implementations at 180 nm technology. Assuming linear scaling,
the RRAM cells follow the scaling trend of SRAM implementations shown in the
below table.

TABLE I: Memories density comparison.

Ref. | Feature Size | Technology | Mb/mm2
---|---|---|---
[23] |  $40\,$nm | CMOS | $0.94$
[24] |  $45\,$nm | CMOS | $0.33$
[3] |  $45\,$nm | FreePDK45 | $0.826$
[25] |  $65\,$nm | CMOS | $0.77$
RRAM Compiler 1 |  $180\,$nm | CMOS | $0.082$
[26]2 |  $180\,$nm | CMOS | $0.067$
[3] |  $0.5\,$µm | SCMOS | $0.005$

1 Memory cell array density. 2 Memory cell size.

##  V Conclusion

To the best of our knowledge, this paper has reported the first open-source
memory compiler specifically for RRAM with supporting RRAM architecture. An
early version of the compiler has been designed and already been made
publically available. Even though, the compiler is a preliminary version, it
is able to generate large sized arrays and verify their layout. The authors of
this document are currently investigating novel suitable circuits to sense the
stored value of a memory cell without altering its resistive state. Next steps
include the research and integration in the compiler of write circuits,
decoders and multiplexers.

Besides providing a powerful tool to academia and industry for boosting memory
design time, the authors believe that the compiler will facilitate the
research of NVMs and RRAMs and their operation.

## Acknowledgement

The authors would like to thank Yihan Pan for providing the 1T1R cell for use
in this memory compiler. This work was in part supported by the Engineering
and Physical Sciences Research Council (EPSRC) Programme under Grant
EP/R024642/1.

## References

  * [1] F. Zahoor _et al._ , “Resistive random access memory (RRAM): an overview of materials, switching mechanism, performance, multilevel cell (MLC) storage, modeling, and applications,” _Nanoscale research letters_ , vol. 15, pp. 1–26, 2020, [Online].
  * [2] A. Chen, “A review of emerging non-volatile memory (NVM) technologies and applications,” _Solid-State Electronics_ , vol. 125, pp. 25–38, 2016, [Online].
  * [3] M. R. Guthaus _et al._ , “Openram: An open-source memory compiler,” in _2016 IEEE/ACM ICCAD_. IEEE, 2016, pp. 1–6, [Online].
  * [4] Y. Xu _et al._ , “A flexible embedded sram ip compiler,” in _2007 IEEE ISCAS_. IEEE, 2007, pp. 3756–3759, [Online].
  * [5] T. Shah _et al._ , “FabMem: A multiported RAM and CAM compiler for superscalar design space exploration.” 2010, [Online].
  * [6] S. Wu _et al._ , “A 65nm embedded low power SRAM compiler,” in _13th IEEE Symposium on DDECS_. IEEE, 2010, pp. 123–124, [Online].
  * [7] R. Goldman _et al._ , “Synopsys’ educational generic memory compiler,” in _10th EWME_. IEEE, 2014, pp. 89–92, [Online].
  * [8] M. Clinton _et al._ , “A 5GHz 7nm L1 cache memory compiler for high-speed computing and mobile applications,” in _2018 IEEE ISSCC_. IEEE, 2018, pp. 200–201, [Online].
  * [9] S. Maheshwari _et al._ , “Hybrid CMOS/memristor circuit design methodology,” 2021. [Online]: http://arxiv.org/abs/2012.02267
  * [10] W. A. Wulf _et al._ , “Hitting the memory wall: Implications of the obvious,” _ACM SIGARCH computer architecture news_ , vol. 23, no. 1, pp. 20–24, 1995, [Online].
  * [11] D. Ielmini _et al._ , “In-memory computing with resistive switching devices,” _Nature Electronics_ , vol. 1, no. 6, pp. 333–343, 2018, [Online].
  * [12] S. Stathopoulos _et al._ , “An electrical characterisation methodology for benchmarking memristive device technologies,” _Scientific reports_ , vol. 9, no. 1, pp. 1–10, 2019, [Online].
  * [13] J. J. Yang _et al._ , “Memristive devices for computing,” _Nature nanotechnology_ , vol. 8, no. 1, pp. 13–24, 2013, [Online].
  * [14] D. B. Strukov _et al._ , “The missing memristor found,” _nature_ , vol. 453, no. 7191, pp. 80–83, 2008, [Online].
  * [15] E. Lee _et al._ , “A ReRAM memory compiler with layout-precise performance evaluation,” in _2019 IEEE SOI-3D-Subthreshold Microelectronics Technology Unified Conference (S3S)_. IEEE, 2019, pp. 1–3, [Online].
  * [16] B. Wicht, _Current Sense Amplifiers for Embedded SRAM in High-Performance System-on-a-Chip Designs: For Embedded SRAM in High-Performance System-on-a-Chip Designs_. Springer Science & Business Media, 2003, vol. 12, [Online].
  * [17] T. N. Blalock _et al._ , “A subnanosecond clamped-bit-line sense amplifier for 1T dynamic RAMs,” in _1991 International Symposium on VLSI-TSA - Proceedings of Technical Papers_. IEEE, 1991, pp. 82–86, [Online].
  * [18] M.-F. Chang _et al._ , “Read circuits for resistive memory (ReRAM) and memristor-based nonvolatile Logics,” in _The 20th ASP-DAC_. IEEE, 2015, pp. 569–574, [Online].
  * [19] _Cadence SKILL Language Reference Product Version ICADVM20.1_ , Cadence Design Systems, San Jose, USA, 2021.
  * [20] _Cadence SKILL Language User Guide Product Version ICADVM20.1_ , Cadence Design Systems, San Jose, USA, 2021.
  * [21] _Calibre® Verification User’s Manual_ , Siemens Mentor, 2020.
  * [22] _Calibre® Interactive™ (New GUI) User’s Manual_ , Siemens Mentor, 2020\.
  * [23] S. Miyano _et al._ , “Highly energy-efficient SRAM with hierarchical bit line charge-sharing method using non-selected bit line charges,” _IEEE JSSC_ , vol. 48, no. 4, pp. 924–931, 2013, [Online].
  * [24] S. O. Toh _et al._ , “Characterization of dynamic SRAM stability in 45 nm CMOS,” _IEEE JSSC_ , vol. 46, no. 11, pp. 2702–2712, 2011, [Online].
  * [25] K. Kushida _et al._ , “A 0.7 v single-supply sram with 0.495 $\mu$m2 cell in 65 nm technology utilizing self-write-back sense amplifier and cascaded bit line scheme,” _IEEE JSSC_ , vol. 44, no. 4, pp. 1192–1198, 2009, [Online].
  * [26] C. A. Kumar _et al._ , “Performance analysis of low power 6T SRAM cell in 180nm and 90nm,” in _2016 2nd International Conference on AEEICB_ , 2016, pp. 351–357, [Online].