# Variability Signatures of a Burst Process in Flaring Gamma-ray Blazars

A. Brill  NASA Postdoctoral Program Fellow, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA

###### Abstract

Blazars exhibit stochastic flux variability across the electromagnetic
spectrum, often exhibiting heavy-tailed flux distributions, commonly modeled
as lognormal. However, Tavecchio et al. (2020) and Adams et al. (2022) found
that the high-energy gamma-ray flux distributions of several of the brightest
flaring Fermi-LAT flat spectrum radio quasars (FSRQs) are well modeled by an
even heavier-tailed distribution, which we show is the inverse gamma
distribution. We propose an autoregressive inverse gamma variability model in
which an inverse gamma flux distribution arises as a consequence of a shot-
noise process. In this model, discrete bursts are individually unresolved and
averaged over within time bins, as in the analysis of Fermi-LAT data.
Stochastic variability on timescales longer than the time bin duration is
modeled using first-order autoregressive structure. The flux distribution
becomes approximately lognormal in the limiting case of many weak bursts. The
fractional variability is predicted to decrease as the time bin duration
increases. Using simulated light curves, we show that the proposed model is
consistent with the typical gamma-ray variability properties of FSRQs and BL
Lac objects. The model parameters can be physically interpreted as the average
burst rate, the burst fluence, and the timescale of long-term stochastic
fluctuations.

††software:  NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007), SciPy
(Virtanen et al., 2020).

##  1 Introduction

Blazars are a class of active galactic nuclei (AGN) which are thought to
possess relativistic jets oriented nearly along our line of sight (Urry &
Padovani, 1995). Blazars make up the vast majority of extragalactic gamma-ray
sources that have been detected in the GeV band by Fermi-LAT (Ajello et al.,
2020).

Blazars are highly variable objects at essentially all wavelengths, with
variability observed in the GeV band at timescales of years down to several
minutes (e.g. Ackermann et al., 2016). Observations of variability can help us
understand the physical processes that give rise to the high-energy emission
from these objects. Characteristic timescales of variability constrain the
apparent size of the emission region by placing upper bounds on the light
crossing time, and are a key input for modeling the physical processes giving
rise to gamma-ray emission.

Blazars can be divided into two source classes, flat spectrum radio quasars
(FSRQs) and BL Lac objects. FSRQs have greater synchrotron luminosity and
lower peak synchrotron frequency than BL Lac objects. Although these two
source classes are often considered to form a continuous “blazar sequence”
(Fossati et al., 1998), they may be better understood as belonging to two
intrinsically distinct AGN populations, of high-power and low-power objects,
respectively (Keenan et al., 2021). In the GeV gamma-ray band, FSRQs are
typically more luminous and more variable than BL Lac objects (Abdo et al.,
2010).

Blazars typically have a power spectral density (PSD) with a power-law shape,
indicating that the emission can be described as a stochastic process, or
random walk (e.g. Abdo et al., 2010). Characteristic variability timescales
may appear as features or spectral breaks in the PSD. Previous works have
found such breaks (when present) at short timescales of $\sim$1 day,
potentially connected to the timescale of successive flare events giving rise
to the high-energy emission, and at long timescales of $\sim$1 year, which may
be related to processes in the accretion disk (e.g. Kataoka et al., 2001; Ryan
et al., 2019; Tarnopolski et al., 2020). However, gamma-ray blazars are also
notable for undergoing intermediate-duration flaring periods which can last
for many days, during which the gamma-ray flux can in some cases increase by
an order of magnitude or more.

One example of the flaring phenomenon comes from the bright FSRQ 3C 279, which
has undergone several such flares (e.g. Ackermann et al., 2016; Adams et al.,
2022). In particular, Adams et al. (2022) analyzed its sub-daily variability
during 10 bright flaring periods lasting between 1 and 11 days, for a total of
54 days of observations. From modeling using exponential profiles, each flare
was resolved into between 1 and 4 separate components, for a total of 24
distinct components that had rise and decay times ranging from days to less
than 1 hour. The fluences of the flare components were found to have a dynamic
range about an order of magnitude narrower than that of their rise and decay
timescales, with a median fluence of $0.85\times 10^{-3}$ erg cm-2.

Blazar variability can be usefully studied in the time domain using
autoregressive processes, a powerful class of time series models well suited
for describing stationary stochastic processes (see Scargle, 1981, and
references therein). In an autoregressive (AR) process of order $p$, the value
$X_{i}$ at a given time step $i$ is a linear combination of the values at the
previous $p$ time steps, to which a stochastic “innovation” term $\epsilon$ is
added. In the related moving average (MA) process of order $q$, a linear
combination of the values of the innovations at the previous $q$ time steps is
used instead. Combining these processes results in a so-called ARMA process of
order $(p,q)$. In this paper, we restrict our attention to AR(1) processes,
that is, the case $p=1$, $q=0$, and we will use the term “autoregressive
process” broadly to refer to the class of ARMA processes in general.

Variability can also be studied through the flux distribution, that is, the
probability density function (PDF) describing the observed flux values.
Blazars are generally found to have heavy-tailed flux distributions, with
excesses at high fluxes compared to a normal distribution. These distributions
are commonly modeled as a lognormal distribution, which can arise through the
multiplicative combination of many independent events (Uttley et al., 2005).
Variability consistent with lognormality has been observed in many studies of
blazars using data in multiple wavebands (e.g. Giebels & Degrange, 2009;
Kushwaha et al., 2017; Valverde et al., 2020; Acciari et al., 2021; Bhatta,
2021).

A commonly observed feature of variable emission from AGN, first discussed in
depth by Uttley & McHardy (2001), is an approximately linear relationship
between the mean flux and standard deviation calculated within subintervals of
the light curve, known as the linear RMS-flux relation. Uttley et al. (2005)
showed that a linear RMS-flux relation can be associated with lognormal
variability, suggesting that the long-term memory needed to produce a linear
RMS-flux relation must be connected with an underlying multiplicative process,
ruling out an additive origin for variability, such as shot noise (Lehto,
1989). Importantly, however, Scargle (2020) showed that an approximately
linear RMS-flux relation arises as a general statistical property of any light
curve that can be modeled as an autoregressive process, with the exact form of
the relation depending on the shape of the innovation. An underlying
multiplicative process is therefore not required to explain the observed
properties of blazar light curves.

In this work, we propose a model in which blazar variability results from an
underlying process of Poisson-distributed bursts. Critically, the shapes of
the individual bursts are taken to be unresolved. The variability then arises
from the statistics of the burst production process, giving rise to an inverse
gamma flux distribution. We derive a model that yields an inverse gamma flux
distribution with AR(1) autoregressive structure. The model has only three
free parameters, representing the average rate of the putative bursts, the
burst fluence, and an autocorrelation timescale representing long-term
fluctuations in the burst rate. These parameters can be interpreted in terms
of physical processes, such as plasmoid-powered flares in a relativistic
magnetic reconnection scenario.

In Section 2, we provide theoretical background on autoregressive processes
and the inverse gamma distribution, and discuss related work. In Section 3, we
motivate the autoregressive inverse gamma model, derive its properties, and
discuss important limiting cases. In Section 4, we use the proposed process to
simulate light curves representative of several source classes. In Section 5,
we discuss some of the model’s limitations, extensions, and implications;
connections to theoretical models; comparisons to previous work; and
applications to blazar classification. Finally, we summarize our conclusions
in Section 6.

##  2 Background and Related Work

An ARMA process of order $(p,q)$ can be written as

$X_{i}=\sum_{j=1}^{p}\phi_{j}X_{i-j}+\sum_{k=1}^{q}\theta_{k}\epsilon_{i-k}+\epsilon_{i},$
(1)

where $X_{i}$ are the values of the time series, $\epsilon_{i}$ are stochastic
innovation terms, and $\phi_{j}$, $\theta_{k}$ are the constant AR and MA
coefficients, respectively. In this paper, we consider the special case of
AR(1) processes with the form

$X_{i}=\phi X_{i-1}+\epsilon_{i}.$ (2)

An AR(1) process in which the innovation terms are assumed to be normally
distributed, such that

$\epsilon\sim\mathcal{N}\left(\mu^{\prime},\sigma^{2}\right),$ (3)

where $\mu^{\prime}\equiv(1-\phi)\mu$ is the mean and $\sigma^{2}$ is the
variance of a Gaussian distribution, is commonly referred to in astrophysics
as an Ornstein-Uhlenbeck or damped random walk process. This process and its
continuous-time analogue have found widespread application in modeling AGN
light curves (e.g. Kelly et al., 2009; Moreno et al., 2019; Burd et al.,
2021). The continuous-time analogue of the general ARMA process, called CARMA,
has been fruitfully deployed as well (Kelly et al., 2014; Ryan et al., 2019).

An AR(1) process is associated with an expected marginal distribution111The
marginal distribution is identical to the expected flux distribution if the
modeled time series is the flux light curve. These distributions will be
different, though closely related, if a different time series is being
considered, such as that of the logarithm of the flux.. If the innovation is
normally distributed according to Eq. 3, the marginal distribution is normally
distributed as well, such that

$X\sim\mathcal{N}\left(\mu,\frac{\sigma^{2}}{1-\phi^{2}}\right).$ (4)

A lognormal flux distribution is usually obtained from an AR(1) model with a
normal innovation by fitting it with the logarithm of the flux, $X=\log F$.
Conveniently, this transformation also prevents the possibility of generating
nonphysical negative fluxes.

One way to interpret the flux distribution is as the result of the combination
of two distinct underlying processes, represented by the autoregressive term
parameterized by $\phi$ and the innovation term parameterized by $\mu$ and
$\sigma$, which may in general have different physical interpretations and
characteristic timescales. The innovation term contributes short-timescale
variability, which is enhanced by the long-timescale variability provided by
the autoregressive term. These are both sources of intrinsic variability,
which are independent of any extrinsic uncertainty contributed by measurement
error.

Recently, Tavecchio et al. (2020) proposed a model of blazar variability based
on a stochastic differential equation (SDE) that exhibits nonlinear dynamics.
Their model is physically motivated by considering a magnetically arrested
accretion disk. The accumulation of magnetic energy is modeled as a
deterministic, equilibrium-reverting process (the autoregressive term) and its
dissipation as a stochastic process (the innovation term). In this model, the
standard deviation of the innovation term depends linearly on the flux, so
that its discrete representation is equivalent to Eq. 2 now with

$\epsilon\sim\mathcal{N}\left(\mu^{\prime},(\sigma X_{i-1})^{2}\right).$ (5)

Assuming that $X>0$, the marginal distribution of this process is given by
(Tavecchio et al., 2020, Eq. A8)

$f(X)=\frac{(\lambda\mu)^{1+\lambda}}{\Gamma(1+\lambda)}\frac{e^{-\lambda\mu/X}}{X^{\lambda+2}},$
(6)

where $\lambda\equiv 2\phi/\sigma^{2}$ and $\Gamma(x)$ is the Gamma
function222In the notation of Tavecchio et al. (2020), $\lambda\equiv
2\theta/\sigma^{2}$.. Tavecchio et al. (2020) showed that this PDF provides a
satisfactory representation of the gamma-ray flux distributions of six bright
FSRQs observed by Fermi-LAT using the light curves analyzed by Meyer et al.
(2019). Subsequently, Adams et al. (2022) demonstrated that this model
provides a better fit than a lognormal PDF to the Fermi-LAT flux distributions
of three bright FSRQs, 3C 279, PKS 1222+216, and Ton 599.

In fact, Eq. 6 has the form of a well-known probability distribution, the
inverse gamma distribution, which has the PDF

$f(X)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\frac{e^{-\beta/X}}{X^{\alpha+1}}.$
(7)

It can be seen that Eq. 6 is equivalent to Eq. 7 with $\alpha=1+\lambda$ and
$\beta=\lambda\mu$.

The inverse gamma distribution and gamma distribution are intimately related.
As the names of these two distributions suggest, if
$X\sim\mathrm{InvGamma}(\alpha,\beta)$, then
$1/X\sim\mathrm{Gamma}(\alpha,\beta)$, and vice-versa333We will also use the
notation $\Gamma(\alpha,\beta)$ with two arguments to denote the Gamma
distribution, and $\Gamma(x)$ with one argument for the Gamma function.. The
PDF of the gamma distribution is given by

$f(X)=\frac{X^{\alpha-1}e^{-\beta X}\beta^{\alpha}}{\Gamma(\alpha)}.$ (8)

The inverse gamma and gamma distributions have positive support, that is,
$X\in(0,\infty)$. For both distributions, $\alpha>0$ is a shape parameter
determining the form of the distribution. Figure 1 shows inverse gamma PDFs
plotted for a range of $\alpha$ values. The inverse gamma distribution is
highly skewed, with a high-flux tail that is heavier than that of a lognormal
distribution and resembles a power law. In the regime $\alpha\gg 1$, the
distribution becomes approximately (but not exactly) lognormal. The scale of
the inverse gamma distribution is determined by $\beta>0$, which for the gamma
distribution is an inverse scale (or rate) parameter.

In a Poisson process, the interarrival time, or time difference between each
event and the next arriving independently, is randomly distributed following
an exponential distribution. The sum of $\alpha$ independent exponential
variables with rate parameter $\beta$ is distributed as
$\mathrm{Gamma}(\alpha,\beta)$. For this reason, the gamma distribution is
closely connected with Poisson processes.

Figure 1: Inverse gamma PDFs for four values of $\alpha$, letting
$\beta=\alpha$. The light thin lines show the lognormal approximation to each
inverse gamma PDF determined using the results derived in Section 3.5. For
$\alpha=100$, the inverse gamma and lognormal curves are indistinguishable on
the scale of the plot.

In this work, we propose an autoregressive inverse gamma model that shares
several commonalities with the model of Tavecchio et al. (2020), including an
inverse gamma flux distribution and autoregressive time structure. However,
the proposed model has a different physical motivation and mathematical basis
deriving from an underlying emission process of discrete unresolved bursts,
rather than from the nonlinear stochastic term of Eq. 5.

These differences lead to several novel predictions. In the proposed model,
variability is controlled primarily by the average burst rate. The model
produces isolated “loner flares” when the average burst rate is low; flaring
light curves reminiscent of highly variable FSRQs and BL Lac objects when it
is moderate; and approximately lognormal variability when it is high. In
addition, high-energy gamma-ray light curves are typically analyzed in time
bins of a set duration. The proposed model predicts that the fractional
variability of a source will decrease proportionately to the square root of
the time bin duration used in the analysis, so long as the bin duration is
longer than the timescale at which individual bursts can be resolved.

##  3 Autoregressive Inverse Gamma Model

###  3.1 Inverse Gamma Flux Distribution

We consider a scenario in which the gamma-ray emission is dominated by
discrete bursts (which might also be called flares or shots) distributed as a
Poisson process. We denote the average arrival rate of the bursts as $r$ and
suppose for simplicity that each burst has the same fluence $\mathcal{F}$. We
consider a light curve binned in time bins of duration $\Delta T$, such that
the observed flux is averaged within each time bin. Therefore, $\Delta T$ is
the minimum timescale at which information about variability in the source can
be determined.

We can roughly associate some characteristic timescale $t_{\mathrm{shape}}$
with the shape of each burst, for example, an exponential rise or decay time.
We make the key assumption that $t_{\mathrm{shape}}\ll\Delta T$, that is, the
bursts are individually unresolved on the timescale of the observations. In
this scenario, all variability driven by the shapes of the individual bursts
is washed out. Instead, the flux distribution derives from the Poisson
statistics of the process producing the bursts.

The average number of bursts occurring in each time bin is given by

$\alpha=r\Delta T.$ (9)

As a first approximation, we assume that each bin has exactly $\alpha$ bursts.
We denote the interarrival time for burst $i$ to occur as
$t_{\mathrm{wait},i}$. Because the observed flux in each time bin is obtained
by averaging over all of the fluxes produced by bursts in that bin, each burst
will contribute a portion of the observed flux given by

$F_{i}=\frac{\mathcal{F}}{t_{\mathrm{wait},i}}.$ (10)

By taking the fluence of burst $i$ to be spread out uniformly on the timescale
$t_{\mathrm{wait},i}$ in its contribution to the observed flux, we have
essentially further assumed that the burst shape is perfectly unresolved. We
discuss the effects of modifying these assumptions in Sections 3.3 and 5.1.

On average, then, each burst will contribute flux $F=\mathcal{F}r$. If the
binned light curve is normalized to units of flux $F_{\mathrm{scale}}$, where
$F_{\mathrm{scale}}$ may be chosen arbitrarily, the average normalized flux of
each burst is

$\beta_{\mathrm{burst}}=\frac{r\mathcal{F}}{F_{\mathrm{scale}}}.$ (11)

From Eq. 10, for burst $i$, the corresponding reciprocal of the flux
$F^{-1}_{i}$ is directly proportional to the interarrival time
$t_{\mathrm{wait},i}$. In a Poisson process, the interarrival time for a
single event is distributed as an exponential distribution. The same is
therefore true for the reciprocal of the flux,
$F^{-1}\sim\mathrm{Exp}(\beta_{\mathrm{burst}})$. Using the relationship
between the exponential and gamma distributions, the reciprocal of the flux
averaged over $\alpha$ bursts in a bin is

$F^{-1}\sim\frac{1}{\alpha}\Gamma(\alpha,\beta_{\mathrm{burst}})=\Gamma(\alpha,\alpha\beta_{\mathrm{burst}})=\Gamma(\alpha,\beta),$
(12)

where $\Gamma(\alpha,\beta)$ is the gamma distribution with shape parameter
$\alpha$ and rate parameter $\beta$, where

$\beta=\alpha\beta_{\mathrm{burst}}=\frac{\alpha
r\mathcal{F}}{F_{\mathrm{scale}}}.$ (13)

Since $F=(F^{-1})^{-1}$, it follows from the relationship between the gamma
and inverse gamma distributions that

$F\sim\mathrm{InvGamma}(\alpha,\beta),$ (14)

where $\mathrm{InvGamma}(\alpha,\beta)$ is the inverse gamma distribution with
shape parameter $\alpha$ and scale parameter $\beta$.

From Eqs. 9 and 13, the parameters of the inverse gamma distribution are
related to the physical characteristics of the bursts as

$\begin{split}r&=\frac{\alpha}{\Delta T}\\\
\mathcal{F}&=\frac{\beta}{\alpha^{2}}F_{\mathrm{scale}}\Delta T.\end{split}$
(15)

###  3.2 Incorporating Autoregression

To model a stationary time series possessing AR(1) autoregressive structure
and an inverse gamma flux distribution, we make use of the existing literature
on conditional linear AR(1) processes (Grunwald et al., 2000), as the simplest
formulation that can produce a non-Gaussian time series with autoregressive
structure. We consider the time series of $F^{-1}$, allowing us to take
advantage of existing models of gamma-distributed AR(1) processes. The flux
light curve, with an inverse gamma flux distribution, can then be obtained by
applying a reciprocal transformation to the time series.

In the standard AR(1) process, a normal innovation produces a normal marginal
distribution, but the gamma distribution does not have this property. A more
complex model is required to produce an AR(1) process with a gamma marginal
distribution. We model the autoregressive structure of the light curve using
the process proposed by Sim (1990), in which the time dependence has the form

$\begin{split}&F_{i}^{-1}=\Gamma(N(F_{i-1}^{-1})+\alpha,\beta),\\\
&N(F^{-1})\sim\mathrm{Pois}(\phi\beta F^{-1}).\end{split}$ (16)

In this approach, known as thinning, rather than multiplying $F_{i-1}^{-1}$ by
a constant $\phi$, it is instead reduced using a stochastic function taking
$\phi$ as a parameter. The process is time-reversible. As shown by Sim (1990),
the autocorrelation function is
$\mathrm{Corr}(F_{i+j}^{-1},F_{i}^{-1})=\phi^{j}$ and the resulting marginal
distribution is

$f(F^{-1})=\Gamma\left(\alpha,(1-\phi)\beta\right).$ (17)

The flux distribution is therefore

$f(F)=\mathrm{InvGamma}(\alpha,(1-\phi)\beta).$ (18)

The autocorrelation parameter $\phi$ can be used to estimate an
autocorrelation timescale $\tau$,

$\phi=e^{-\Delta T/\tau},$ (19)

or as a function of $\phi$,

$\tau=\frac{\Delta T}{\ln{(1/\phi)}}.$ (20)

Several other models for a gamma AR(1) process have been proposed, possessing
different statistical properties. In the GAR(1) model (Gaver & Lewis, 1980;
Lawrance, 1982; Walker, 2000), a marginal gamma distribution is obtained by
considering an innovation constructed to represent a shot-noise process with
an exponential distribution of shot amplitudes. The process is not time-
reversible and exhibits peaks followed by geometrically decaying runs of
values. While we do not consider this model well suited to describe fluxes
averaged over time bins containing multiple unresolved bursts, it may have
useful applications in other astrophysical contexts in which fast rising and
exponentially decaying flares are individually resolved.

In the BGAR(1) process (Lewis et al., 1989), the innovation has a gamma
distribution, and a gamma marginal distribution is obtained by replacing the
constant autoregression coefficient $\phi$ with a random variable following a
beta distribution with expected value $\phi$. The BGAR(1) process is time-
reversible. However, the marginal distribution resulting from this process is
independent of $\phi$, which is inconsistent with the dependence on $\phi$ in
the marginal distribution that might be expected by considering the limiting
Gaussian case (Eq. 4). Here, the autoregressive structure could be thought of
as making it likely that nearby time steps sample similar values of the
innovation PDF, which is difficult to interpret physically.

###  3.3 Accounting for the Discretization of Bursts in Time Bins

As we have described the model so far, each time bin contains exactly the
average number of bursts $\alpha$. More realistically, the number of bursts
would vary from time bin to time bin following a Poisson distribution with
parameter $\alpha$. This scenario has two main implications. First, a time bin
could contain no bursts, so it would have no associated flux. To allow this,
we can describe the inverse flux using a mixed distribution, with finite
probability mass $p_{0}=e^{-\alpha}$ at $F^{-1}=0$ and continuous probability
density integrating to $p_{\Gamma}=1-e^{-\alpha}$ for $F^{-1}>0$.

Figure 2: Best-fit inverse gamma flux distributions to Poisson-weighted
mixtures of inverse gamma distributions, using the approximation method
described in Sec. 3.3.

Second, for bins that do have bursts, because the number of bursts is
different in each bin, the continuous part of the inverse-flux distribution is
no longer a single gamma distribution with shape parameter $\alpha$, but a
weighted mixture of gamma distributions with shape parameters k = 1, 2, 3,
etc. While the flux distribution resulting from this mixture is not exactly
inverse gamma in general, as shown in Figure 2, it can be approximated well by
a single inverse gamma distribution.

The properties of the resulting process are derived in Appendix A. The
distribution of $F^{-1}>0$ is given by a gamma distribution with shape
parameter

$\alpha_{\mathrm{obs}}=\frac{A(\alpha)}{1-e^{-\alpha}}\alpha,$ (21)

and rate parameter

$\beta_{\mathrm{obs}}=A(\alpha)(1-\phi)\beta,$ (22)

where

$A(\alpha)=\left(\frac{\alpha}{(1-e^{-\alpha})^{2}}\sum_{l=1}^{\infty}\frac{\alpha^{l}e^{-\alpha}}{l!}\frac{1}{l}\right)^{-1}.$
(23)

Figure 3: Top: Modification factor $A(\alpha)$ of the inverse gamma parameters
to account for a Poisson distribution of bursts, as a function of the
intrinsic average burst rate $\alpha$. Bottom: The apparent burst rate
$\alpha_{\mathrm{obs}}$ as a function of $\alpha$.

In Figure 3, $A(\alpha)$ and $\alpha_{\mathrm{obs}}$ are plotted as a function
of $\alpha$. As shown in the figure, $A(\alpha)\lesssim 1$ for all $\alpha$,
reaching a minimum of $A(\alpha)\approx 0.73$ at $\alpha\approx 2.98$, and
approaching 1 in the limits $\alpha\to 0$ and $\alpha\to\infty$. As a result,
$\alpha_{\mathrm{obs}}$ approaches 1 as $\alpha\to 0$ and approaches $\alpha$
as $\alpha\to\infty$.

For time bins containing bursts, the behavior of the time series is described
by Eq. 16 with a simple change of parameters. Since $\alpha_{\mathrm{obs}}$
and $\beta_{\mathrm{obs}}$ are the parameters of the marginal distribution,
the process parameters are $\alpha_{\mathrm{obs}}$ and
$\beta_{\mathrm{obs}}/(1-\phi)=A(\alpha)\beta$, giving

$\begin{split}&F_{i}^{-1}=\Gamma(N(F_{i-1}^{-1})+\alpha_{\mathrm{obs}},A(\alpha)\beta),\\\
&N(F^{-1})\sim\mathrm{Pois}(\phi A(\alpha)\beta F^{-1}).\end{split}$ (24)

To accommodate bins with $F^{-1}=0$ while preserving the AR(1) autocorrelation
structure, the modified process can be described as a Markov process with two
states, $F^{-1}=0$ for bins without bursts and
$F^{-1}\sim\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}})$ for bins with
bursts. The state transition probabilities $s_{i}\in\\{0,\Gamma\\}\to
s_{i+1}\in\\{0,\Gamma\\}$ are given by

$\begin{split}p_{00}&=1-p(\phi,\alpha)(1-e^{-\alpha})\\\
p_{0\Gamma}&=p(\phi,\alpha)(1-e^{-\alpha})\\\ p_{\Gamma
0}&=p(\phi,\alpha)e^{-\alpha}\\\
p_{\Gamma\Gamma}&=1-p(\phi,\alpha)e^{-\alpha},\end{split}$ (25)

where

$p(\phi,\alpha)=\frac{1-\phi}{1+\frac{\phi}{\alpha}A^{-1}(\alpha)(1-e^{-\alpha})}.$
(26)

###  3.4 “Loner Flares” at Low Burst Rates

When the burst rate is very low, most time bins contain no bursts, and the
light curve consists of isolated flares separated by intervals of no emission.
The model parameters $\alpha$, $\beta$, and $\phi$ can be extracted from the
properties of the isolated flares by considering the limiting case where
$\alpha\ll 1$. Specifically, to first order in $\alpha$, the parameters of the
observed distribution become

$\displaystyle\alpha_{\mathrm{obs}}$ $\displaystyle\approx
1+\frac{\alpha}{4},$ (27a) $\displaystyle\beta_{\mathrm{obs}}$
$\displaystyle\approx\left(1-\frac{\alpha}{4}\right)(1-\phi)\beta,$ (27b)
$\displaystyle p_{\Gamma}$ $\displaystyle\approx\alpha,$ (27c) $\displaystyle
p(\phi,\alpha)$
$\displaystyle\approx\frac{1-\phi}{1+\phi}\left(1+\frac{\phi}{1+\phi}\frac{\alpha}{4}\right).$
(27d)

The expected fraction of time bins that contain flaring gamma-ray emission
equals $p_{\Gamma}$ and so provides a direct estimate of $\alpha$ from Eq.
27c. The average length of flares in units of $\Delta T$ is given by

$N_{\mathrm{avg}}=\frac{1}{p_{\Gamma
0}}\approx\frac{1+\phi}{1-\phi}\left(1+\frac{\alpha}{1+\phi}\right).$ (28)

In particular, if $\alpha$ is very small, $\phi$ can be roughly estimated as

$\phi\approx\frac{N_{\mathrm{avg}}-1}{N_{\mathrm{avg}}+1}.$ (29)

A measurement of $\beta_{\mathrm{obs}}$, thereby giving an estimate of $\beta$
from Eq. 27b, can then be obtained by fitting an inverse gamma distribution to
the observed flaring flux distribution.

###  3.5 Lognormal Behavior at High Burst Rates

When the burst rate is very high, the law of large numbers applies: every time
bin contains a number of bursts close to the expected value, and the dynamic
range of variability is reduced. The inverse flux distribution approaches a
normal distribution, resulting in a flux distribution that is approximately
lognormal. Specifically, if $\alpha$ is large, the parameters of the observed
distribution become $\alpha_{\mathrm{obs}}\approx\alpha$,
$\beta_{\mathrm{obs}}\approx(1-\phi)\beta$, and $p_{\Gamma}\approx 1$, and

$\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}})\to\mathcal{N}\left(\frac{\alpha_{\mathrm{obs}}}{\beta_{\mathrm{obs}}},\frac{\alpha_{\mathrm{obs}}}{\beta_{\mathrm{obs}}^{2}}\right),$
(30)

where $\mathcal{N}(\mu,\sigma^{2})$ is the normal distribution with mean $\mu$
and variance $\sigma^{2}$. Then, to a level of $s$ standard deviations, the
flux $F$ is constrained such that

$\left|\frac{1}{F}-\frac{\alpha_{\mathrm{obs}}}{\beta_{\mathrm{obs}}}\right|\lesssim\frac{s\sqrt{\alpha_{\mathrm{obs}}}}{\beta_{\mathrm{obs}}},$
(31)

from which we obtain

$\left|\frac{\beta_{\mathrm{obs}}}{\alpha_{\mathrm{obs}}}\frac{1}{F}-1\right|\lesssim\frac{s}{\sqrt{\alpha}}.$
(32)

Letting $\Delta\equiv\beta_{\mathrm{obs}}/(\alpha_{\mathrm{obs}}F)-1$, then
within a significance level of $\sim 5$ standard deviations,
$\left|\Delta\right|\lesssim 1$ when $\alpha\gtrsim 25$.

For a random variable $X$ with PDF $f(X)$, the reciprocal variable $Y=1/X$ has
the PDF $g(Y)=Y^{-2}f(1/Y)$. Since the inverse flux is normally distributed,
the flux $F$ has the PDF

$\begin{split}F&\sim\frac{1}{F^{2}}\frac{1}{(\sqrt{\alpha_{\mathrm{obs}}}/\beta_{\mathrm{obs}})\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{(1/F)-(\alpha_{\mathrm{obs}}/\beta_{\mathrm{obs}})}{\sqrt{\alpha_{\mathrm{obs}}}/\beta_{\mathrm{obs}}}\right)^{2}}\\\
&=(1+\Delta)\frac{1}{F}\sqrt{\frac{\alpha_{\mathrm{obs}}}{2\pi}}e^{-\frac{1}{2}\alpha_{\mathrm{obs}}\Delta^{2}}\\\
&\approx\frac{1}{F}\sqrt{\frac{\alpha_{\mathrm{obs}}}{2\pi}}e^{-\frac{1}{2}\alpha_{\mathrm{obs}}\ln^{2}{(1+\Delta)}}\\\
&=\frac{1}{F}\sqrt{\frac{\alpha_{\mathrm{obs}}}{2\pi}}e^{-\frac{1}{2}\alpha_{\mathrm{obs}}(\ln{F}-\ln{(\beta_{\mathrm{obs}}/\alpha_{\mathrm{obs}}}))^{2}}\\\
&=\mathrm{Lognormal}(\mu,\sigma^{2}),\end{split}$ (33)

where

$\begin{split}\mu=\ln{\frac{\beta_{\mathrm{obs}}}{\alpha_{\mathrm{obs}}}}&\approx\ln{\frac{(1-\phi)\beta}{\alpha}},\\\
\sigma^{2}=\frac{1}{\alpha_{\mathrm{obs}}}&\approx\frac{1}{\alpha}.\end{split}$
(34)

To first order, the flux distribution can be approximated as a lognormal
distribution, the parameters of which can be interpreted in terms of the
underlying burst process. Using Eq. 15, the burst parameters can be estimated
from the lognormal fit as

$\begin{split}r&\approx\left(\sigma^{2}\Delta T\right)^{-1}\\\
\mathcal{F}&\approx\frac{1}{1-\phi}\sigma^{2}e^{\mu}F_{\mathrm{scale}}\Delta
T.\end{split}$ (35)

###  3.6 Fractional Variability

The shape of the flux distribution depends on the shape parameter $\alpha$.
From Eq. 9, $\alpha$ is the product of two factors: the burst rate $r$, which
is an intrinsic physical parameter, and the time bin duration $\Delta T$,
which is an extrinsic choice made in the analysis. This fact enables us to
predict how the fractional variability of a given source should scale when
analyzed using time bins of different durations.

The mean of an inverse gamma distribution is

$\mu=\frac{\beta}{\alpha-1},$ (36)

and its variance is

$\sigma^{2}=\frac{\beta^{2}}{(\alpha-1)^{2}(\alpha-2)}.$ (37)

The mean is defined only for $\alpha>1$ and the variance for $\alpha>2$. Then,
for $\alpha>2$, the expected fractional variability is

$\begin{split}F_{\mathrm{var}}&=\frac{\sigma}{\mu}=(\alpha_{\mathrm{obs}}-2)^{-1/2}\approx(\alpha-2)^{-1/2}\\\
&=r^{-1/2}\left(1-\frac{2}{r\Delta T}\right)^{-1/2}\Delta
T^{-1/2}.\end{split}$ (38)

The fractional variability of a given source should scale with the time
binning as $F_{\mathrm{var}}\propto\Delta T^{-1/2}$, at least in the regime
$\Delta T\gtrsim\max(2/r,t_{\mathrm{shape}})$.

##  4 Representative Light Curves

Figure 4: Top Left: Simulated light curve meant to resemble a “loner flare”
blazar. Top Right: Simulated light curve meant to resemble an FSRQ, with
parameters chosen based on gamma-ray observations of 3C 279. Bottom Left:
Simulated light curve meant to resemble a “flaring BL Lac” object, such as 1ES
1215+303. Bottom Right: Simulated light curve meant to resemble a “constant BL
Lac” object. Figure 5: Flux distributions corresponding to the four simulated
light curves. The expected PDF is shown as an orange line in each plot. Figure
6: RMS-flux relations corresponding to the four simulated light curves. A
linear least-squares fit is shown with a dotted black line.

To illustrate the process, we show simulated 12-year weekly light curves
representing four different source types in Figure 4. The light curves were
generated following the formulas given in Eqs. 24 and 25. Each light curve
represents only one possible realization of the corresponding stochastic
process.

To generate a representative light curve for an FSRQ, we estimated $r$ and
$\mathcal{F}$ from the observations reported by Adams et al. (2022) of flares
of 3C 279, giving $r\approx 24/54~{}\mathrm{bursts}~{}\mathrm{day}^{-1}\approx
5\times 10^{-6}~{}\mathrm{s}^{-1}$ and $\mathcal{F}\approx 0.85\times
10^{-3}~{}\mathrm{erg}~{}\mathrm{cm}^{-2}$. We set a fairly high
autocorrelation parameter of $\phi=0.95$, consistent with the absence of a
resolvable spectral break in the gamma-ray PSD of 3C 279 as determined both by
CARMA modeling (Ryan et al., 2019) and by the power spectral response method
(Goyal et al., 2021). With these parameters, $\alpha_{\mathrm{obs}}\approx
2.3$, similar to the best-fit values of $\lambda+1=1.7$ and 2.0 reported
respectively by Tavecchio et al. (2020) and Adams et al. (2022).

The resulting simulated light curve is shown in the top right of Figure 4. The
pattern of variability, showing several flaring periods separated by intervals
of relatively quiescent activity, indeed appears to resemble an actual light
curve of an FSRQ like 3C 279. The overall normalization of the energy flux
during the flares and quiescent periods is also similar to that seen in the
actual light curve of that source.

Figure 4, top left, illustrates a source exhibiting even more extreme
variability, with isolated flares separated by intervals of complete
quiescence. For this simulation, the burst rate was decreased by a factor of
50 and the burst fluence increased by the same. This light curve resembles
that of a so-called “loner flare” blazar characterized by the detection of a
single high-flux flaring event in the long-term gamma-ray light curve (Wang et
al., 2020).

The autoregressive inverse gamma model may also be able to reproduce the less
extreme variability associated with BL Lac objects. For example, an extensive
analysis of the multiwavelength variability of the flaring high-synchrotron-
peaked BL Lac object 1ES 1215+303 was performed by Valverde et al. (2020), who
found that its flux distribution was described better by a lognormal
distribution than by a normal distribution. They obtained a lognormal
distribution with a best-fit value of $\sigma=0.43\pm 0.02$, which, using Eq.
34, can be reproduced by an inverse gamma distribution with $\alpha\approx
5.4$.

With this motivation, we generated simulated light curves with parameters
chosen to represent two “BL Lac” objects, a “flaring” source with
$\alpha_{\mathrm{obs}}\approx 5$ intended to be reminiscent of a source like
1ES 1215+303, and a “constant” source with $\alpha_{\mathrm{obs}}\approx 60$.
These simulated light curves are shown in Figure 4, bottom left and right,
respectively. Again, the simulated light curves appear to resemble real light
curves.

The flux distributions corresponding to the simulated light curves shown in
Figure 4 are plotted in Figure 5, along with the expected PDF for each.

We examined the RMS-flux relation for each of the simulated light curves. The
light curves were binned in intervals of 25 data points and the mean and
standard deviation was calculated within each interval. The standard deviation
for each interval is plotted against the mean in Figure 6. An approximately
linear RMS-flux relation can be observed in each case. The “constant BL Lac”
light curve, with the least skewed flux distribution, exhibits the RMS-flux
relation with the most scatter.

##  5 Discussion

###  5.1 Applicability and Extensions of the Model

The autoregressive inverse gamma model for blazar variability presented in
this work is based on several simplifying assumptions which should be
considered when interpreting the model and the resulting flux distribution. An
in-depth exploration of the effects of modifying these assumptions is beyond
the scope of the simple statistical model presented in this work, and may be
better addressed in the context of detailed physical simulations.

First, we have assumed that all bursts have exactly the same fluence. More
realistically, the fluence may vary from burst to burst. An illustration of
how varying the fluence can affect the flux distribution is shown in Figure 7.
Distributions are shown for $\alpha=1$, 2, 5 and 60, similar to the simulated
sources considered in Section 4, scaled such that $\beta=\alpha$. We
numerically simulated a modification of the basic scenario described in Sec.
3.1, without autoregression, for variable burst fluence by multiplying the
mean fluence $\mathcal{F}$ by a random scale factor for each burst. The scale
factors were drawn from a truncated Gaussian distribution with $\mu=1$ and
$\sigma=0.275$. The distribution was truncated at $-3\sigma$ to avoid negative
fluctuations and at $+3\sigma$ to preserve symmetry. We chose a value of
$\sigma=0.275$ to yield a fluence distribution varying by one order of
magnitude.

Figure 7: Simulated flux distributions incorporating variations in burst
fluence of up to an order of magnitude. The corresponding unmodified inverse
gamma distributions are shown with orange lines.

Allowing the fluence to vary produces an excess at the low-flux end of the
flux distribution, causing it to cut off less sharply, while preserving the
heavy high-flux tail. This occurs because the fluence fluctuations are
dominated by the burst arrival-time fluctuations except at the lowest fluxes,
at least for small $\alpha$. For the much narrower distribution with
$\alpha=60$, the entire distribution is noticeably shifted, while retaining an
approximately lognormal shape.

An additional low-flux component could also arise in the flux distribution for
other reasons. Such a component might come about, for example, through the
flux contributed by a subdominant population of weak bursts with low fluence.
Also, if the bursts have relatively long rise or decay times, “leakage” of the
flux from one bin to the next could result in a similar effect. The effect of
such a component would be most apparent as an excess at the low end of the
flux distribution, in addition to potentially adding a small flux component to
the time bins that are considered here to have no flux. However, an additional
low-flux component may be difficult to identify in data, as it must be
disentangled from potential artifactual contributions, such as statistical
fluctuations or imperfectly modeled gamma-ray flux coming from nearby sources
in the field of view.

In addition, as shown in Figure 2, the exact Poisson-weighted mixture
distribution has in general heavier tails than the best-fit inverse gamma
distribution. For this reason, a pure inverse gamma flux distribution may
somewhat underestimate the true expected variability.

We have also assumed that the burst arrival times follow a Poisson process,
but in reality there may be physical constraints on how often bursts can
occur. For example, if a lower bound exists on the time that must elapse
between two bursts, that would place an upper bound on the imputed flux of a
single burst as well as on the number of bursts that can occur in a time bin,
inducing a cut-off at the high-flux end of the flux distribution.

The high-flux end of the flux distribution is critical for distinguishing the
inverse gamma distribution from other proposed distributions, such as the
lognormal. However, high-flux time bins are rare. The longer the light curve,
the more effectively the high-flux part of the flux distribution can be
studied. Some previous studies that compared normal and lognormal
distributions have excluded flaring intervals from the flux distribution so as
not to overly favor the heavier-tailed lognormal distribution (e.g. Giebels &
Degrange, 2009; Valverde et al., 2020). However, when only considering heavy-
tailed distributions, these intervals should not be excluded.

The derivation of the statistical model presented here requires that
individual bursts are individually unresolved. For this reason, the model is
only applicable at time scales longer than the putative burst shape timescale
$t_{\mathrm{shape}}$. For very bright sources, it may be possible to estimate
$t_{\mathrm{shape}}$ from an analysis of short-timescale variability during
high-amplitude flares, but for weaker sources, this quantity may be more
difficult to measure. At timescales shorter than $t_{\mathrm{shape}}$, the
observed variability would depend on the characteristics of the bursts as well
as on the statistics of their stochastic production process.

For example, in the study done by Adams et al. (2022), a time binning of 1 day
was used when studying the flux distribution of 3C 279, which is shorter than
the longest timescale of $\sim$100 hr found by fitting exponential profiles to
individual flare components. In this case, care must be taken when
interpreting the inverse gamma fit parameters using the model presented in
this work.

In this paper, we exclusively considered AR(1) processes. Some gamma-ray
blazars may exhibit possible year-scale (quasi-)periodicities (e.g. Peñil et
al., 2020; Rueda et al., 2022) or trends (e.g. Valverde et al., 2020). Future
work could investigate the possibility of extending the model by adding
deterministic trend or periodic terms, higher-order autoregressive terms, or
moving average terms. Sources with possible year-scale quasi-periodic
oscillations could be modeled using a second-order autoregressive process.
Alternatively, one could study nonlinear time-series structures in which the
innovation depends on the flux, such as in the SDE proposed by Tavecchio et
al. (2020) or the autoregressive gamma process (ARG) of Gourieroux & Jasiak
(2006).

###  5.2 Connection to Theoretical Models

The simulations shown in Section 4 suggest that a single stochastic process
can describe both flares and quiescent periods in light curves, consistent
with the previous studies that have modeled blazar light curves using
autoregressive processes (e.g. Ryan et al., 2019). The process proposed in
this work has three free parameters, connected to the average burst rate $r$,
the typical burst fluence $\mathcal{F}$, and the characteristic
autocorrelation timescale $\tau$. These quantities can be interpreted in terms
of the physical processes producing the gamma-ray emission.

One theoretical scenario that may give rise to a burst process of gamma-ray
emission is magnetic reconnection, in which compact magnetized structures
called plasmoids are produced that give rise to fast (minutes - days) gamma-
ray flares (Giannios, 2013). Magnetic reconnection models based on particle-
in-cell simulations predict that large, slow-moving plasmoids and much smaller
but relativistic ones both produce flares with comparable fluence (Petropoulou
et al., 2016). This finding suggests a motivation for the assumption applied
in this work that all bursts have approximately the same fluence.

Simulations of plasmoid growth and mergers have found that in a single
reconnection event, the bolometric light curve is dominated by flares from a
few large plasmoids, along with a superposition of very fast variability
contributed by smaller plasmoids (Petropoulou et al., 2018; Christie et al.,
2019). Fast variability is especially evident in the high-energy gamma-ray
light curve (Acciari et al., 2020) in comparison to other wave bands, such as
optical (Zhang et al., 2022). Simulations of magnetic reconnection have been
shown to yield flux levels and variability patterns compatible with gamma-ray
observations of flaring blazars, particularly FSRQs (Meyer et al., 2021).

The autocorrelation timescale is most likely associated with longer-timescale
variability. It may be directly connected to the same underlying physics as
the burst process, such as very long duration flares produced by large
plasmoids in a magnetic reconnection scenario (Giannios, 2013), or it may be
indirectly related and associated with a different physical process. In
particular, the variability of the jet emission may originate in processes in
the accretion disk. In this case, $\tau$ may reflect one of several timescales
associated with dynamical, thermal, or viscous processes in the disk (Czerny,
2006). Alternatively, quasi-periodic variability with a characteristic
timescale of approximately 1 to 100 days could result from a “striped jet”
characterized by magnetic fields of alternating polarity (Giannios & Uzdensky,
2019). The variability of the light curve may be driven by different physical
processes at different time scales. If so, $\tau$ may change when the light
curve is examined using different choices of $\Delta T$.

###  5.3 Comparison to Previous Work

The model presented in this paper shares several important similarities with
the SDE proposed by Tavecchio et al. (2020), along with a number of key
differences. The model of Tavecchio et al. (2020) is motivated by the
consideration of processes in a magnetically arrested accretion disk. In the
model proposed in this work, the form of the flux distribution instead derives
from the statistics of a burst process giving rise to the gamma-ray emission,
although the long-term variations captured by the autoregressive structure may
be considered to be related to processes in the accretion disk.

Both models predict AR(1) time structure and flux distributions with the form
of an inverse gamma distribution. However, the models exhibit different
behavior when the shape parameter of the distribution becomes small. In the
autoregressive inverse gamma model, as $\alpha\to 0$,
$\alpha_{\mathrm{obs}}\to 1$ and
$\beta_{\mathrm{obs}}\to(1-\phi)\beta_{\mathrm{burst}}$. Only a few time bins
have nonzero flux, resulting in a light curve characteristic of a loner flare
blazar. In the model of Tavecchio et al. (2020), in the limiting case of
$\lambda\to 0$ indicating the predominance of the stochastic term, the
resulting flux distribution approaches a pure power-law distribution,
$f(X)\propto X^{-2}$.

In addition, the model proposed by Tavecchio et al. (2020) is represented as
an SDE. The process is fundamentally continuous, although discrete
approximations are numerically useful for generating simulations and fitting
the model to data. The quality of the discrete approximation would be expected
to increase for data that are increasingly finely binned in time (up to the
effective time resolution of the instrument). On the other hand, the model
proposed in this work is inherently discrete in character, as it is derived
from a stochastic process of individual bursts. At short timescales, it would
be expected to break down as the time structure of the bursts themselves
begins to dominate the light curve. Mathematically, this discreteness comes
about because the innovation is not Gaussian, which would be required to
obtain a continuous Markov process (Gillespie, 1996).

In the model proposed by Tavecchio et al. (2020), which is only supposed to be
a simplified representation of the real dynamics, the stochastic term consists
of Gaussian fluctuations multiplied by the normalized flux. As such,
stochastic fluctuations can drive the flux below zero, especially when the
flux is large. The requirement that the flux be greater than zero is required
as an implicit additional constraint for deriving an inverse gamma flux
distribution. In the model proposed in this work, the flux is inherently
nonnegative, since the probability distributions have positive support.

The autoregressive inverse gamma model predicts a heavy-tailed flux
distribution and approximately linear RMS-flux relation, similar to lognormal
variability models based on multiplicative processes (Uttley et al., 2005).
However, it explains these phenomena with an additive burst process modulated
by long-term stochastic variability, without invoking multiplicative
processes. In this way, the proposed model is more similar in its physical
interpretation to shot-noise models that explain the light curve as resulting
from the additive superposition of many overlapping bursts (e.g. Lehto, 1989;
Tanihata et al., 2001). However, there is a significant difference: because
the bursts are unresolved, the flux distribution and power spectrum are
determined by the statistics of the burst arrival process and the long-term
stochastic variability, rather than by the shapes and amplitudes of the bursts
themselves.

Duda & Bhatta (2021) studied gamma-ray blazar variability using the log-stable
family of probability distributions, which generalizes the lognormal
distribution. Most of the blazars in their sample displayed variability
consistent with lognormality; intriguingly, however, a subset of the sources
had extremely heavy-tailed flux distributions leading to infinite variance.
This finding of infinite variance is characteristic of the inverse gamma
distribution with $\alpha_{\mathrm{obs}}<2$. Interestingly, the sample of
sources with infinite variance examined by Duda & Bhatta (2021) contained both
BL Lac objects and FSRQs.

Another model class that produces a heavy-tailed marginal distribution is that
of generalized autoregressive conditional heteroskedasticity (GARCH) models,
commonly used in finance for modeling time series of price fluctuations
(Engle, 1982; Bollerslev, 1986; Francq & Zakoian, 2019). GARCH models describe
white-noise processes in which the variance of the innovation changes over
time, following an ARMA process. The time-varying variance is called the
volatility.

The autoregressive inverse gamma model proposed in this work differs
significantly from GARCH models, as it has a fixed innovation distribution.
The model of Tavecchio et al. (2020), which incorporates volatility in its
stochastic term, is related to GARCH models. However, it differs from a
standard first-order GARCH model by also having a mean-reverting
autoregressive term.

The GARCH(1, 1) process has an inverse gamma marginal distribution of
volatility, giving the resulting time series a Student’s t marginal
distribution (Nelson, 1990). Golan & Gerig (2013) proposed an alternative
model for autocorrelated inverse gamma volatility in price fluctuations, in
which exogenous events influencing trading activity play a similar role to the
burst process considered in this work. The model proposed in this work goes
beyond that of Golan & Gerig (2013) by allowing the number of bursts (or
events) to fluctuate from time bin to time bin, including the possibility of
no bursts (Sec. 3.3); enforcing an exact AR(1) autocorrelation structure; and
having separate parameters modeling the underlying burst rate and
autocorrelation timescale.

###  5.4 Blazar Classification

In the autoregressive inverse gamma model, the extreme variability typical of
bright flaring FSRQs, such as that studied by Tavecchio et al. (2020) and
Adams et al. (2022), is associated with a small burst rate $r$. Relating high-
energy gamma-ray variability to physically interpretable quantities may help
characterize the differences observed between FSRQs and BL Lac objects,
increasing our understanding of the relationship between these blazar classes.
For example, Uemura et al. (2020) identified differences between the optical
light curves of FSRQs and BL Lac objects by modeling them using Ornstein-
Uhlenbeck processes. If FSRQs and BL Lac objects make up distinct classes
rather than a continuous sequence, we might expect the burst parameters of
these objects to belong to distinct clusters, with the opposite being true
otherwise.

In some cases, blazar classifications are uncertain or disputed. For example,
Padovani et al. (2019) have argued that the blazar TXS 0506+056, which has
been associated with the detection of high-energy neutrinos, is an object
intrinsically of the FSRQ type masquerading as a BL Lac object. Gamma-ray
variability modeling may provide evidence to help us better understand these
sources.

##  6 Summary and Conclusions

Models of gamma-ray variability offer a critical tool for understanding the
physical processes producing high-energy emission in blazars, potentially
helping us characterize the relationship between FSRQs and BL Lac objects. In
this work, we have proposed an autoregressive inverse gamma model of blazar
variability that provides a simple framework for interpreting gamma-ray blazar
variability on multiple timescales. In the proposed model, an inverse gamma
flux distribution derives from an emission process in which gamma rays are
produced in discrete bursts arriving as a Poisson process. Importantly, the
bursts are taken to be individually unresolved within the time bins used in
the analysis; only the average flux of each time bin is observed. Long-
timescale stochastic variations are modeled by incorporating first-order
autoregressive, or AR(1), structure.

The autoregressive inverse gamma model has three free parameters, representing
the average burst rate, the burst fluence, and the autocorrelation timescale.
These parameters can be interpreted in terms of physical quantities, such as
by associating the bursts with plasmoid-powered flares in a magnetic
reconnection scenario, for example. Furthermore, in the proposed model, flares
and quiescent periods lasting from days to months are caused by random
fluctuations in the arrival rate of the bursts. This intermediate-timescale
activity naturally emerges from the interaction between two underlying
physical processes, the short-timescale burst process and the long-timescale
stochastic variations. Flaring and quiescent emission therefore have the same
origin.

The autoregressive inverse gamma model yields simulated gamma-ray light curves
consistent with multiple blazar source classes. The variability
characteristics are controlled primarily by the burst rate. In particular, the
use of parameters estimated from flare observations of the FSRQ 3C 279 yields
a simulated light curve that demonstrates the variability behavior
characteristic of that source, with bright flares separated by longer
relatively quiescent periods, along with a similar overall flux normalization.
Decreasing the burst rate generates a light curve consisting of isolated
“loner flares”, while increasing it produces light curves more closely
resembling those of BL Lac objects. In this case, flaring activity is less
prominent, and the flux distribution is approximately lognormal.

The proposed model is based on an emission process of discrete bursts,
resulting in novel predictions that differentiate it from previous work,
including the model proposed by Tavecchio et al. (2020) that also features an
inverse gamma flux distribution and autoregressive time structure. In
particular, the fractional variability is predicted to decrease as a larger
time binning is used, scaling as $F_{\mathrm{var}}\propto\Delta T^{-1/2}$.

Because the model has AR(1) structure, the PSD has the form of a power law
with a low-frequency break. A high-frequency break may also occur, depending
on the characteristics of the unresolved burst process.

Thanks to Alberto Domínguez, Reshmi Mukherjee, Jeremy Perkins, Jeff Scargle,
and Janeth Valverde for helpful discussions, comments, and suggestions on this
work. A.B. is supported by the NASA Postdoctoral Program at Goddard Space
Flight Center, administered by USRA and ORAU.

##  Appendix A Derivation of an AR(1) Gamma Process with Poisson-distributed
Burst Counts

In this section, we re-derive the gamma process proposed by Sim (1990) to
incorporate a more realistic model of burst counts. To do so, we replace the
shape parameter $\alpha=r\Delta T$ that represents the average number of
bursts in each bin with a Poisson distribution with parameter $\alpha$.

Since $\alpha$ is Poisson distributed, we have

$\alpha_{\mathrm{bin}}=\begin{cases}0,&p=e^{-\alpha}\\\
l,&p=\frac{\alpha^{l}e^{-\alpha}}{l!},\end{cases}$ (A1)

for natural numbers $l>0$. A randomly chosen bin will with probability
$e^{-\alpha}$ contain no bursts, and therefore generate no flux. It follows
that time bins with observed flux will always have $\alpha_{\mathrm{bin}}>0$,
which is also a requirement for the process of Sim (1990) to be well-defined.

We therefore wish to construct an AR(1) process that switches between a zero
and non-zero state (we will show later that in fact $F^{-1}=0$ for bins with
no bursts). To do so, we take inspiration from the Binary AR(1) process
(McKenzie, 1985). Writing $F^{-1}\equiv X$ in the remainder of this section
for brevity, we consider a Markov process with two states, $X=0$ and
$X\sim\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}})$, with transition
probabilities for $X_{i}\to X_{i+1}$ notated as

$\begin{cases}p_{00},&0\to 0\\\
p_{0\Gamma},&0\to\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}})\\\
p_{\Gamma 0},&\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}})\to 0\\\
p_{\Gamma\Gamma},&\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}})\to\Gamma(\alpha_{\mathrm{obs}},\beta_{\mathrm{obs}}).\end{cases}$
(A2)

To ensure a self-consistent process, we have the relations

$\displaystyle p_{00}+p_{0\Gamma}$ $\displaystyle=1$ (A3a) $\displaystyle
p_{\Gamma 0}+p_{\Gamma\Gamma}$ $\displaystyle=1$ (A3b) $\displaystyle
p_{0}+p_{\Gamma}$ $\displaystyle=1,$ (A3c)

where $p_{0}=e^{-\alpha}$ and $p_{\Gamma}=1-e^{-\alpha}$, as well as

$\displaystyle p_{0}p_{00}+p_{\Gamma}p_{\Gamma 0}$ $\displaystyle=p_{0}$ (A4a)
$\displaystyle p_{0}p_{0\Gamma}+p_{\Gamma}p_{\Gamma\Gamma}$
$\displaystyle=p_{\Gamma}.$ (A4b)

######  Remark.

Let an isolated flare be defined as a contiguous sequence of bins with
$X\sim\Gamma$ surrounded by bins with $X=0$ on both ends. The expected length
of such a sequence is given by the expected number of time bins needed for one
$\Gamma\to 0$ transition to occur, starting from an initial time bin with
$X\sim\Gamma$. The average isolated flare length in units of number of time
bins is therefore

$L=\frac{1}{p_{\Gamma 0}}.$ (A5)

######  Remark.

Let the effective duration of a time bin with $X\sim\Gamma$ be defined as the
total duration encompassing that time bin and any subsequent time bins with
$X=0$, until the next bin with $X\sim\Gamma$ is reached. The expected length
of a sequence of bins with $X=0$ is $1/p_{0\Gamma}$. Making use of this fact
and Eqs. A4b and A3c, the expected effective duration is

$N_{\mathrm{eff}}=\Delta T\left(1+\frac{p_{\Gamma
0}}{p_{0\Gamma}}\right)=\Delta
T\left(1+\frac{p_{0}}{p_{\Gamma}}\right)=\frac{\Delta
T}{p_{\Gamma}}=\frac{\Delta T}{1-e^{-\alpha}}.$ (A6)

We now focus on the case $\alpha_{\mathrm{bin}}>0$, that is, the process that
takes place when the transition $\Gamma\to\Gamma$ occurs, and will return
later to finish considering the remaining cases. Renormalizing the probability
distribution of $\alpha_{\mathrm{bin}}$ for $\alpha_{\mathrm{bin}}>0$ yields

$\alpha_{\mathrm{bin,obs}}=\begin{cases}l,&p=\frac{1}{1-e^{-\alpha}}\frac{\alpha^{l}e^{-\alpha}}{l!}.\end{cases}$
(A7)

For a time bin with $\alpha_{\mathrm{bin,obs}}=l$, the corresponding observed
rate is

$\beta_{\mathrm{bin,obs}}=\frac{\Delta
T}{N_{\mathrm{eff}}}\beta_{\mathrm{burst}}l=\frac{1-e^{-\alpha}}{\alpha}\beta
l,$ (A8)

where two modifications have been incorporated. First, the flux is averaged
over $l$ bursts rather than $\alpha$ bursts. Second, the average flux is
reduced by a “leakage” factor of $\Delta T/N_{\mathrm{eff}}$ because the
fluence produced in a time bin is supposed to be averaged over the burst
interarrival time, which is the entire effective duration $N_{\mathrm{eff}}$,
but the flux is only collected within that bin itself, of duration $\Delta T$.

The process mean $\mathbb{E}(X)$ and variance $\mathrm{Var}(X)$ can be
calculated as follows. From the law of total expectation, we have

$\begin{split}\mathbb{E}(X)&=\mathbb{E}(X_{i+1})\\\
&=\mathbb{E}_{X_{i}}\left[\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{1}{1-e^{-\alpha}}\frac{\alpha^{l}e^{-\alpha}}{l!}\mathbb{E}_{X_{i+1}}[\Gamma(X_{i+1};k+l,\frac{1-e^{-\alpha}}{\alpha}\beta
l)]\right]\\\
&=\frac{1}{(1-e^{-\alpha})^{2}}\frac{\alpha}{\beta}\mathbb{E}_{X_{i}}\left[\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{\alpha^{l}e^{-\alpha}}{l!}\left(\frac{k}{l}+1\right)\right]\\\
&=\frac{1}{(1-e^{-\alpha})^{2}}\frac{\alpha}{\beta}\mathbb{E}_{X_{i}}\left[\sum_{l=1}^{\infty}\frac{\alpha^{l}e^{-\alpha}}{l!}\left(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
X_{i}+1\right)\right]\\\
&=\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta}\mathbb{E}_{X_{i}}\left[\phi\frac{1-e^{-\alpha}}{\alpha}\beta
X_{i}+1\right]\\\
&=\phi\mathbb{E}(X)+\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta}\\\
&=\frac{1}{1-\phi}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta},\end{split}$
(A9)

where we have made use of the fact that for a stationary process,
$\mathbb{E}(X_{i+1})=\mathbb{E}(X_{i})=\mathbb{E}(X)$.

From the law of total variance, we have

$\begin{split}\mathrm{Var}(X)&=\mathbb{E}(\mathrm{Var}(X_{i+1}|X_{i}))+\mathrm{Var}(\mathbb{E}(X_{i+1}|X_{i}))\\\
&=\mathbb{E}(\mathrm{Var}(X_{i+1}|X_{i}))+\mathrm{Var}\left(\phi
X_{i}+\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta}\right)\\\
&=\mathbb{E}(\mathrm{Var}(X_{i+1}|X_{i}))+\phi^{2}\mathrm{Var}(X)\\\
&=\frac{1}{1-\phi^{2}}\mathbb{E}(\mathrm{Var}(X_{i+1}|X_{i})).\end{split}$
(A10)

For a mixture of distributions $f_{i}$ with means $\mu_{i}$, variances
$\sigma^{2}_{i}$, and weights $p_{i}$, the variance of the mixture $f=\sum
p_{i}f_{i}$ is given by

$\mathrm{Var}(f)=\sum_{i}p_{i}\left(\sigma_{i}^{2}+\mu_{i}^{2}\right)-\left(\sum_{i}p_{i}\mu_{i}\right)^{2}.$
(A11)

Using Eq. A11, we can therefore write

$\begin{split}\begin{aligned} \mathrm{Var}(X_{i+1}|X_{i})=\>\\\
\end{aligned}&\overbrace{\begin{aligned}
\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{1}{1-e^{-\alpha}}\frac{\alpha^{l}e^{-\alpha}}{l!}\biggl{[}&\mathrm{Var}(\Gamma(X_{i+1};k+l,\frac{1-e^{-\alpha}}{\alpha}\beta
l))\\\ &+\mathbb{E}(\Gamma(X_{i+1};k+l,\frac{1-e^{-\alpha}}{\alpha}\beta
l))^{2}\biggr{]}\end{aligned}}^{(A)}\\\
&-\underbrace{\left(\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{1}{1-e^{-\alpha}}\frac{\alpha^{l}e^{-\alpha}}{l!}\mathbb{E}(\Gamma(X_{i+1};k+l,\frac{1-e^{-\alpha}}{\alpha}\beta
l))\right)^{2}}_{(B)},\end{split}$ (A12)

where

$\begin{split}(A)&=\frac{1}{1-e^{-\alpha}}\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{\alpha^{l}e^{-\alpha}}{l!}\frac{(k+l)+(k+l)^{2}}{\left(\frac{1-e^{-\alpha}}{\alpha}\beta
l\right)^{2}}\\\
&=\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}\frac{\alpha}{(1-e^{-\alpha})^{2}}\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{\alpha^{l}e^{-\alpha}}{l!}\frac{k^{2}+2kl+l^{2}+k+l}{l^{2}}\\\
&=\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}\frac{\alpha}{(1-e^{-\alpha})^{2}}\sum_{l=1}^{\infty}\frac{\alpha^{l}e^{-\alpha}}{l!}\left[\left(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
X_{i}\right)^{2}+2\phi\frac{1-e^{-\alpha}}{\alpha}\beta
X_{i}+1+\frac{2\phi\frac{1-e^{-\alpha}}{\alpha}\beta X_{i}+1}{l}\right]\\\
&=(\phi X_{i})^{2}+\frac{2\phi\alpha\beta
X_{i}}{1-e^{-\alpha}}+\left(\frac{\alpha}{1-e^{-\alpha}}\right)^{2}+\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}\left(2\phi\frac{1-e^{-\alpha}}{\alpha}\beta
X_{i}+1\right)A^{-1}(\alpha),\end{split}$ (A13)

with

$A(\alpha)=\left(\frac{\alpha}{(1-e^{-\alpha})^{2}}\sum_{l=1}^{\infty}\frac{\alpha^{l}e^{-\alpha}}{l!}\frac{1}{l}\right)^{-1},$
(A14)

and

$\begin{split}(B)&=\frac{1}{(1-e^{-\alpha})^{4}}\frac{\alpha^{2}}{\beta^{2}}\left[\sum_{k=0}^{\infty}\sum_{l=1}^{\infty}\frac{(\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i})^{k}e^{-\phi\frac{1-e^{-\alpha}}{\alpha}\beta
lX_{i}}}{k!}\frac{1}{1-e^{-\alpha}}\frac{\alpha^{l}e^{-\alpha}}{l!}\left(\frac{k}{l}+1\right)\right]^{2}\\\
&=\left(\phi
X_{i}+\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta}\right)^{2}\mathrm{(From~{}Eq.~{}\ref{eq:mu})}\\\
&=(\phi X_{i})^{2}+\frac{2\phi\alpha\beta
X_{i}}{1-e^{-\alpha}}+\left(\frac{\alpha}{1-e^{-\alpha}}\right)^{2}.\end{split}$
(A15)

Combining Eqs. A10, A12, A13, and A15 and making use of Eq. A9 in the second
step,

$\begin{split}\mathrm{Var}(X)&=\frac{1}{1-\phi^{2}}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}A^{-1}(\alpha)\mathbb{E}\left[2\phi\frac{1-e^{-\alpha}}{\alpha}\beta
X_{i}+1\right]\\\
&=\frac{1}{1-\phi^{2}}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}A^{-1}(\alpha)\left(2\phi\frac{1-e^{-\alpha}}{\alpha}\beta\frac{1}{1-\phi}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta}+1\right)\\\
&=\frac{1}{1-\phi^{2}}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}A^{-1}(\alpha)\left(\frac{2\phi}{1-\phi}+1\right)\\\
&=\frac{1}{(1-\phi)^{2}}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta^{2}}A^{-1}(\alpha).\end{split}$
(A16)

The resulting mixture of gamma distributions does not have the exact form of a
gamma distribution, but it can be reasonably well approximated as one and
approaches an exact gamma distribution in the limits of $\alpha\to 0$ and
$\alpha\to\infty$. The best-fit parameters of the mixture to a gamma
distribution can be estimated using the method of moments,

$\alpha_{\mathrm{obs}}=\frac{\mathbb{E}(X)^{2}}{\mathrm{Var}(X)}=\frac{A(\alpha)}{1-e^{-\alpha}}\alpha,$
(A17)

and

$\beta_{\mathrm{obs}}=\frac{\mathbb{E}(X)}{\mathrm{Var}(X)}=A(\alpha)(1-\phi)\beta.$
(A18)

We can now validate our assumption that $F^{-1}=0$ for bins with no bursts.
The mean of the average process $X_{\mathrm{avg}}\sim\Gamma(\alpha,\beta)$ is

$\mathbb{E}(X_{\mathrm{avg}})=\frac{1}{1-\phi}\frac{\alpha}{\beta},$ (A19)

which should be equal to the mean of the combined process,

$\begin{split}\mathbb{E}(X_{\mathrm{avg}})&=e^{-\alpha}\mathbb{E}(X_{0})+(1-e^{-\alpha})\mathbb{E}(X_{\mathrm{obs}})\\\
\frac{1}{1-\phi}\frac{\alpha}{\beta}&=e^{-\alpha}\mathbb{E}(X_{0})+(1-e^{-\alpha})\frac{1}{1-\phi}\frac{1}{1-e^{-\alpha}}\frac{\alpha}{\beta},\end{split}$
(A20)

from which we obtain $\mathbb{E}(X_{0})=0$. Since $F^{-1}\geq 0$, it follows
that $F^{-1}=0$ for all bins with no bursts. The flux for bins with no bursts
is therefore undefined.

We now return to the full process, considering both bins with and without
bursts. To impose on the process AR(1) correlation structure with
autocorrelation constant $\phi$, we use the condition

$\phi=\frac{\mathrm{Cov}(X_{i+1},X_{i})}{\mathrm{Var}(X)}=\frac{\mathbb{E}(X_{i+1}X_{i})-\mathbb{E}(X)^{2}}{\mathrm{Var}(X)}.$
(A21)

From the law of total expectation,

$\begin{split}\mathbb{E}(X_{i+1}X_{i})&=p_{0}p_{00}\cdot
0^{2}+(p_{0}p_{0\Gamma}+p_{\Gamma}p_{\Gamma 0})\cdot
0\cdot\mathbb{E}(X_{\Gamma})+p_{\Gamma}p_{\Gamma\Gamma}\mathbb{E}(X_{\Gamma,i},X_{\Gamma,i+1})\\\
&=p_{\Gamma}p_{\Gamma\Gamma}\mathbb{E}(X_{\Gamma,i},X_{\Gamma,i+1})\\\
&=p_{\Gamma}p_{\Gamma\Gamma}\left[\phi\mathrm{Var}(X_{\Gamma})+\mathbb{E}(X_{\Gamma})^{2}\right]\\\
&=p_{\Gamma\Gamma}\frac{1}{(1-\phi)^{2}}\frac{\alpha^{2}}{\beta^{2}}\left(\frac{\phi}{\alpha
A(\alpha)}+\frac{1}{1-e^{-\alpha}}\right)\\\ \end{split}$ (A22)

Using Eq. A11,

$\begin{split}\mathrm{Var}(X)&=p_{0}\left(\mathrm{Var}(X_{0})+\mathbb{E}(X_{0})^{2}\right)+p_{\Gamma}\left(\mathrm{Var}(X_{\Gamma})+\mathbb{E}(X_{\Gamma})^{2}\right)-\left(p_{0}\mathbb{E}(X_{0})+p_{\Gamma}\mathbb{E}(X_{\Gamma})\right)^{2}\\\
&=p_{\Gamma}\left[\mathrm{Var}(X_{\Gamma})+\mathbb{E}(X_{\Gamma})^{2}-p_{\Gamma}\mathbb{E}(X_{\Gamma})^{2}\right]\\\
&=\frac{1}{(1-\phi)^{2}}\frac{\alpha^{2}}{\beta^{2}}\left(\frac{1}{\alpha
A(\alpha)}+\frac{e^{-\alpha}}{1-e^{-\alpha}}\right).\end{split}$ (A23)

Combining Eqs. A21, A22, A19, and A23 and simplifying yields

$p_{\Gamma\Gamma}=1-\frac{(1-\phi)e^{-\alpha}}{1+\frac{\phi}{\alpha}A^{-1}(\alpha)(1-e^{-\alpha})}=1-p(\phi,\alpha)e^{-\alpha},$
(A24)

where

$p(\phi,\alpha)=\frac{1-\phi}{1+\frac{\phi}{\alpha}A^{-1}(\alpha)(1-e^{-\alpha})}.$
(A25)

From Eqs. A3a and A3b, and Eq. A4a or A4b, we obtain the full set of
transition probabilities,

$\begin{split}p_{00}&=1-p(\phi,\alpha)(1-e^{-\alpha})\\\
p_{0\Gamma}&=p(\phi,\alpha)(1-e^{-\alpha})\\\ p_{\Gamma
0}&=p(\phi,\alpha)e^{-\alpha}\\\
p_{\Gamma\Gamma}&=1-p(\phi,\alpha)e^{-\alpha}.\end{split}$ (A26)

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 722, 520, doi: 10.1088/0004-637X/722/1/520
  * Acciari et al. (2020) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2020, ApJS, 248, 29, doi: 10.3847/1538-4365/ab89b5
  * Acciari et al. (2021) —. 2021, MNRAS, 504, 1427, doi: 10.1093/mnras/staa3727
  * Ackermann et al. (2016) Ackermann, M., Anantua, R., Asano, K., et al. 2016, ApJ, 824, L20, doi: 10.3847/2041-8205/824/2/L20
  * Adams et al. (2022) Adams, C. B., Batshoun, J., Benbow, W., et al. 2022, ApJ, 924, 95, doi: 10.3847/1538-4357/ac32bd
  * Ajello et al. (2020) Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e
  * Bhatta (2021) Bhatta, G. 2021, ApJ, 923, 7, doi: 10.3847/1538-4357/ac2819
  * Bollerslev (1986) Bollerslev, T. 1986, Journal of Econometrics, 31, 307, doi: 10.1016/0304-4076(86)90063-1
  * Burd et al. (2021) Burd, P. R., Kohlhepp, L., Wagner, S. M., et al. 2021, A&A, 645, A62, doi: 10.1051/0004-6361/202039097
  * Christie et al. (2019) Christie, I. M., Petropoulou, M., Sironi, L., & Giannios, D. 2019, MNRAS, 482, 65, doi: 10.1093/mnras/sty2636
  * Czerny (2006) Czerny, B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 360, AGN Variability from X-Rays to Radio Waves, ed. C. M. Gaskell, I. M. McHardy, B. M. Peterson, & S. G. Sergeev, 265
  * Duda & Bhatta (2021) Duda, J., & Bhatta, G. 2021, MNRAS, 508, 1446, doi: 10.1093/mnras/stab2574
  * Engle (1982) Engle, R. F. 1982, Econometrica, 50, 987, doi: 10.2307/1912773
  * Fossati et al. (1998) Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433, doi: 10.1046/j.1365-8711.1998.01828.x
  * Francq & Zakoian (2019) Francq, C., & Zakoian, J.-M. 2019, GARCH Models: Structure, Statistical Inference and Financial Applications (John Wiley & Sons), doi: 10.1002/9781119313472
  * Gaver & Lewis (1980) Gaver, D. P., & Lewis, P. 1980, Advances in Applied Probability, 12, 727, doi: 10.2307/1426429
  * Giannios (2013) Giannios, D. 2013, MNRAS, 431, 355, doi: 10.1093/mnras/stt167
  * Giannios & Uzdensky (2019) Giannios, D., & Uzdensky, D. A. 2019, MNRAS, 484, 1378, doi: 10.1093/mnras/stz082
  * Giebels & Degrange (2009) Giebels, B., & Degrange, B. 2009, A&A, 503, 797, doi: 10.1051/0004-6361/200912303
  * Gillespie (1996) Gillespie, D. T. 1996, American Journal of Physics, 64, 225, doi: 10.1119/1.18210
  * Golan & Gerig (2013) Golan, R., & Gerig, A. 2013, arXiv e-prints, arXiv:1306.4975.  https://arxiv.org/abs/1306.4975
  * Gourieroux & Jasiak (2006) Gourieroux, C., & Jasiak, J. 2006, Journal of Forecasting, 25, 129, doi: 10.1002/for.978
  * Goyal et al. (2021) Goyal, A., Soida, M., Stawarz, L., et al. 2021, arXiv e-prints, arXiv:2112.04194.  https://arxiv.org/abs/2112.04194
  * Grunwald et al. (2000) Grunwald, G. K., Hyndman, R. J., Tedesco, L., & Tweedie, R. L. 2000, Australian & New Zealand Journal of Statistics, 42, 479, doi: 10.1111/1467-842X.00143
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kataoka et al. (2001) Kataoka, J., Takahashi, T., Wagner, S. J., et al. 2001, ApJ, 560, 659, doi: 10.1086/322442
  * Keenan et al. (2021) Keenan, M., Meyer, E. T., Georganopoulos, M., Reddy, K., & French, O. J. 2021, MNRAS, 505, 4726, doi: 10.1093/mnras/stab1182
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kelly et al. (2014) Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., & Uttley, P. 2014, ApJ, 788, 33, doi: 10.1088/0004-637X/788/1/33
  * Kushwaha et al. (2017) Kushwaha, P., Sinha, A., Misra, R., Singh, K. P., & de Gouveia Dal Pino, E. M. 2017, ApJ, 849, 138, doi: 10.3847/1538-4357/aa8ef5
  * Lawrance (1982) Lawrance, A. 1982, Scandinavian Journal of Statistics, 234
  * Lehto (1989) Lehto, H. J. 1989, in ESA Special Publication, Vol. 1, Two Topics in X-Ray Astronomy, Volume 1: X Ray Binaries. Volume 2: AGN and the X Ray Background, ed. J. Hunt & B. Battrick, 499
  * Lewis et al. (1989) Lewis, P. A., McKenzie, E., & Hugus, D. K. 1989, Stochastic Models, 5, 1, doi: 10.1080/15326348908807096
  * McKenzie (1985) McKenzie, E. 1985, Journal of the American Water Resources Association, 21, 645, doi: 10.1111/j.1752-1688.1985.tb05379.x
  * Meyer et al. (2021) Meyer, M., Petropoulou, M., & Christie, I. M. 2021, ApJ, 912, 40, doi: 10.3847/1538-4357/abedab
  * Meyer et al. (2019) Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39, doi: 10.3847/1538-4357/ab1651
  * Moreno et al. (2019) Moreno, J., Vogeley, M. S., Richards, G. T., & Yu, W. 2019, PASP, 131, 063001, doi: 10.1088/1538-3873/ab1597
  * Nelson (1990) Nelson, D. B. 1990, Journal of Econometrics, 45, 7, doi: 10.1016/0304-4076(90)90092-8
  * Padovani et al. (2019) Padovani, P., Oikonomou, F., Petropoulou, M., Giommi, P., & Resconi, E. 2019, MNRAS, 484, L104, doi: 10.1093/mnrasl/slz011
  * Peñil et al. (2020) Peñil, P., Domínguez, A., Buson, S., et al. 2020, ApJ, 896, 134, doi: 10.3847/1538-4357/ab910d
  * Petropoulou et al. (2018) Petropoulou, M., Christie, I. M., Sironi, L., & Giannios, D. 2018, MNRAS, 475, 3797, doi: 10.1093/mnras/sty033
  * Petropoulou et al. (2016) Petropoulou, M., Giannios, D., & Sironi, L. 2016, MNRAS, 462, 3325, doi: 10.1093/mnras/stw1832
  * Rueda et al. (2022) Rueda, H., Glicenstein, J.-F., & Brun, F. 2022, ApJ, 934, 6, doi: 10.3847/1538-4357/ac771c
  * Ryan et al. (2019) Ryan, J. L., Siemiginowska, A., Sobolewska, M. A., & Grindlay, J. 2019, ApJ, 885, 12, doi: 10.3847/1538-4357/ab426a
  * Scargle (1981) Scargle, J. D. 1981, ApJS, 45, 1, doi: 10.1086/190706
  * Scargle (2020) —. 2020, ApJ, 895, 90, doi: 10.3847/1538-4357/ab8d38
  * Sim (1990) Sim, C. H. 1990, Journal of Applied Probability, 27, 325, doi: 10.2307/3214651
  * Tanihata et al. (2001) Tanihata, C., Urry, C. M., Takahashi, T., et al. 2001, ApJ, 563, 569, doi: 10.1086/323960
  * Tarnopolski et al. (2020) Tarnopolski, M., Żywucka, N., Marchenko, V., & Pascual-Granado, J. 2020, ApJS, 250, 1, doi: 10.3847/1538-4365/aba2c7
  * Tavecchio et al. (2020) Tavecchio, F., Bonnoli, G., & Galanti, G. 2020, MNRAS, 497, 1294, doi: 10.1093/mnras/staa2055
  * Uemura et al. (2020) Uemura, M., Abe, T., Yamada, Y., & Ikeda, S. 2020, PASJ, 72, 74, doi: 10.1093/pasj/psaa063
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Uttley & McHardy (2001) Uttley, P., & McHardy, I. M. 2001, MNRAS, 323, L26, doi: 10.1046/j.1365-8711.2001.04496.x
  * Uttley et al. (2005) Uttley, P., McHardy, I. M., & Vaughan, S. 2005, MNRAS, 359, 345, doi: 10.1111/j.1365-2966.2005.08886.x
  * Valverde et al. (2020) Valverde, J., Horan, D., Bernard, D., et al. 2020, ApJ, 891, 170, doi: 10.3847/1538-4357/ab765d
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walker (2000) Walker, S. G. 2000, Scandinavian Journal of Statistics, 27, 575, doi: 10.1111/1467-9469.00208
  * Wang et al. (2020) Wang, G., Wang, Z., Chen, L., Zhou, J., & Xing, Y. 2020, PASJ, 72, 9, doi: 10.1093/pasj/psz128
  * Zhang et al. (2022) Zhang, H., Li, X., Giannios, D., et al. 2022, ApJ, 924, 90, doi: 10.3847/1538-4357/ac3669