# Sachdev-Ye-Kitaev model: Non-self-averaging properties of the energy
spectrum

Richard Berkovits  Department of Physics, Jack and Pearl Resnick Institute,
Bar-Ilan University, Ramat-Gan 52900, Israel

###### Abstract

The short time (large energy) behavior of the Sachdev-Ye-Kitaev model (SYK) is
one of the main motivation to the growing interest garnered by this model.
True chaotic behavior sets in at the Thouless time, which can be extracted
from the energy spectrum. In order to do so, it is necessary to unfold the
spectrum, i.e., to filter out global tendencies. Using a simple ensemble
average for unfolding results in a parametically low estimation of the
Thouless energy. By examining the behavior of the spectrum as the distribution
of the matrix elements is changed into a log-normal distribution it is shown
that the sample to sample level spacing variance determines this estimation of
the Thouless energy. Using the singular value decomposition method, SVD, which
filters out these sample to sample fluctuations, the Thouless energy becomes
parametrically much larger, essentially of order of the band width. It is
shown that the SYK model in non-self-averaging even in the thermodynamic limit
which must be taken into account in considering its short time properties.

##  I Introduction

The interplay between disorder and interactions in quantum systems has been a
central theme in condensed-matter physics for the last half century. Recently
the Sachdev-Ye-Kitaev (SYK) model [1, 2], has garnered much interest in the
fields of quantum gravity and quantum field theory [3, 4]. A key feature of
the model is that it follows random matrix theory (RMT) behavior, as is
manifested in the chaotic behavior of its energy spectrum.

The SYK model first appeared in the context of spin liquids [1] and then in
string theory [5] and quantum gravity [4]. The SYK model is known to be
chaotic [6, 7, 8], showing a Wigner-like behavior of the energy spectra. Once
the SYK model is perturbed by a single-body random term which mimics diagonal
disorder in the Anderson model [9, 10, 11, 12, 13, 14, 15], or several SYK
dots are coupled by single-body random terms [16, 17, 18, 19, 20], a
transition from metallic (chaotic) to insulating behavior occurs which leads
to a Wigner to Poisson crossover of the statistical properties of the spectra.

While studying nuclear and condensed-matter systems it became clear that many
physical systems exhibits universal behavior at long times (short energy
scales) [21, 22, 23, 24, 25, 26]. Nevertheless, universality may break at
shorter times (large energy scales) for which a particle had no time to sample
the entire phase space of the system and its behavior depends on local non-
universal features [27]. Thus, in the context of disordered metals the scale
for which the metallic spectrum deviates from the universal behavior is known
as the Thouless energy, $E_{Th}=\hbar D/{\tilde{L}}^{2}=g\Delta$ ($D$ is the
diffusion constant, ${\tilde{L}}$, is the linear dimension, $g$ is the
dimensionless conductance, $\Delta$ the average level spacing) and the
Thouless time $t_{Th}=\hbar/E_{Th}={\tilde{L}}^{2}/D$.

The question whether an analogue of the Thouless energy manifests itself in
the SYK model has surfaced in Ref. 6 where García-García and Verbaarschot have
studied (among other things) the variance of the number of levels as function
of the size of an energy window $E$, denoted by $\langle\delta^{2}n(E)\rangle$
(where $\langle\ldots\rangle$ represents an ensemble average and $n(E)$ is the
number of levels within $E$). A departure from the RMT behavior is apparent
above a certain value of $E$, which quite naturally was identified with the
Thouless energy. Moreover, at larger energy windows, the number variance
adopts a quadratic form $\langle\delta^{2}n(E)\rangle=a+b\langle
n(E)\rangle^{2}$. Evidence for the Thouless energy has also been seen for
other measures such as off diagonal expectation values [28]. In Ref. 19 the
origin of the Thouless energy for the SYK was identified as the relaxation of
modes prevailing at shorter scale. For disordered metals these modes are know
as the diffusion modes. For SYK it was suggested that similar modes can be
constructed, where the number of such modes is connected to the number of of
independent interaction terms in the SYK model which is much smaller than the
size of the Hilbert space. This leads to an energy scale which determines the
Thouless energy. Thus, the energy scale for the departure of the SYK model
level number variance from the RMT prediction is determined by the scale of
the sample to sample fluctuations of the ensemble.

For a single disordered realization of the SYK model all the terms have just a
single variance scale, which is well behaved (box or Gaussian) and therefore
the origin of the additional energy scale determining the Thouless energy is
not obvious. As noted recently by Jia and Verbaarschot [30] this deviation
from RMT behavior has its root in large scale sample to sample fluctuations of
the spectrum. They attribute these fluctuations to the relatively small number
of independent random variables contributing to the SYK Hamiltonian. When
these small number of long-wavelength sample to sample fluctuations are
parameterized by terms of Q-Hermite orthogonal polynomials and removed from
the spectrum of a particular realization a pure RMT behavior is retained up to
a very large energy scale.

The influence could be quantified [30] by estimating the energy scale for
which the small number of independent random variables will change the number
variance. Using the notation for the complex SYK (CSYK) half-filled model
defined in the next section, where $L$ is the number of sites, the size of the
Hilbert space is $\binom{L}{L/2}\sim 2^{L}$. The number of independent
variables is $\binom{2L}{4}\sim L^{4}$ leading to a variance of $L^{-2}$ in
any observable. Thus one expects the number of levels $n$ to deviate
significantly $O(1)$ from RMT predictions on a scale of $n\sim L^{2}$. A
similar result emerges from the calculation in Ref. 19 where the coefficient
$b$ in the number variance was estimated as $b\sim L^{-4}$ thus becoming
significant at $n\sim L^{2}$.

Although the deviation from RMT is shifted to larger $n$ as the system size
$L$ increases, the proportion of levels following RMT predictions out of the
the total number of states goes to zero as $L^{2}/2^{L}$. On the other hand,
one would expect that after filtering out long wavelength sample to sample
fluctuations the RMT behavior will be followed for a finite portion of the
spectrum. Thus, one expects the Thouless energy to crucially depend on whether
one simply averages over an ensemble or takes into account the single sample
adjustments. This is a hallmark of a non-self-averaging system [31]. The
number variance for the ensemble average is different than the number variance
adjusted to a particular sample.

Hence, the energy scale for which the spectrum departs from the RMT
predictions crucially depends on the the unfolding, i.e., the method by which
the averaged over the density of states is performed. Estimating the local
density of states by a simple ensemble average will give a different value
than an unfolding method that is able to take into account sample specific
global behavior of the spectrum. In recent studies it has been shown that [32,
33, 34, 35, 36, 37, 38], these sample specific long ranged features of the
spectrum can be identified by the singular value decomposition (SVD)
procedure. As detailed below, similar in a sense to Fourier transform, SVD
actually reconstructs the energy spectrum of each realization by a sum over a
series of SVD amplitudes multiplied by the corresponding SVD mode. Unlike the
Fourier transform, the amplitudes of the SVD are identical for all
realizations while the modes are realization specific. Generally, plotting the
amplitudes from large to small (known as a Scree plot) shows that the largest
amplitudes (usually $O(1)$ modes) are orders of magnitude larger than the
rest, while the following modes amplitudes obey a power law. The largest
amplitudes modes depict the very long wave length behavior of the spectra,
while the modes whose amplitudes follow a power law capture the shorter range
properties. Thus SVD is a very natural method to examine the SYK spectrum
behavior, and uncover the realization specific universal properties of the
spectrum.

One may conclude that there are two possible energy scales for the deviation
of the spectrum of the SYK model from RMT predictions. The first is the energy
scale corresponding to the deviation from the ensemble average, which in self
averaging systems such as disordered metals is the Thouless energy. The second
energy scale corresponds to the deviation from RMT when the sample to sample
long-wavelength fluctuations are removed from the spectrum. This sample
specific parameterization of these long range behavior can be done using
Q-Hermite orthogonal polynomials [30] or by SVD [32, 33, 34, 35, 36, 37, 38],
or probably by other method. The energy scale of these fluctuations is another
relevant energy scale which for self-averaging systems such as disordered
metals is equivalent to the Thouless energy obtained from the ensemble average
[36]. For a non-self-averaging system such as the SYK model, these two energy
scales are not equal, and for clarity we will retain the notation of Thouless
energy, $E_{Th}$, for the case where the spectrum is unfolded by a local
average over all realizations, while the realization adjusted Thouless energy
$E_{Th^{*}}$ is obtained using a realization adapted unfolding method.

Here we will show that one can tweak the behavior of the CSYK model to a more
non-self-averaging behavior by changing the distribution of the off-diagonal
to a log-normal distribution. Thus it is possible to enhance the sample to
sample fluctuations and study its influence on $E_{Th}$ and $E_{Th^{*}}$. Such
a wide distribution was not previously considered for the SYK model and should
help to clarify the divergence of $E_{Th^{*}}$ from $E_{Th}$.

The paper is arranged as follows: CSYK is defined in Sec. II. Corroborating
the expected behavior for short range energy scales is performed in Sec. III.
The universal statistics 4-fold symmetry as function of the system size is
observed. In Sec. IV long range energy scales are probed by the number
variance. The spectrum is unfolded by using the local ensemble averaged level
spacing. RMT predictions hold only up to $E_{Th}$, which becomes smaller as
the log-normal distribution acquires a thicker tail towards larger values.
Above $E_{Th}$ the variance increases quadratically as function of the number
of levels in the energy window. In Sec. V we switch to the SVD analysis. This
analysis reveals that each realization has a distinct long-range correlated
level spacing structure. Thus, one should adapt the unfolding for each
realization. After a realization adapted unfolding the number variation
follows the RMT prediction for much larger energy scales, i.e., $E_{Th^{*}}\gg
E_{Th}$. Actually, $E_{Th}$ could be estimated from the sample to sample
variance in the level spacing. In Sec. VI these results are discussed in the
limit of large CSYK systems, showing the in the thermodynamic limit the SYK is
non-self-averaging.

##  II Complex SYK model

Here we use the the complex spinless fermions version of the SYK model given
by the following Hamiltonian:

$\displaystyle\hat{H}=\sum_{i>j>k>l}^{L}V_{i,j,k,l}\hat{c}_{i}^{\dagger}\hat{c}_{j}^{\dagger}\hat{c}_{k}\hat{c}_{l},$
(1)

the couplings $V_{i,j,k,l}$ are complex numbers, where the real and imaginary
components are independently drawn from an identical distribution. We study
here two different distribution: The first is a box distribution between
$-L^{-3/2}/2\ldots L^{-3/2}/2$, where L is the number of sites. The second
distribution is the log-normal distribution

$\displaystyle P(V)=(A/|V|)e^{\frac{-\ln^{2}(|V|/V_{\rm typ})}{2p\ln(V_{\rm
typ}^{-1})}},$ (2)

with $V_{\rm typ}\sim K^{-\gamma/2}$, $K$ is the size of the Hilbert space,
$A$ a normalization and $\gamma$, $p$ are parameters. For simplicity here we
shall set $p=0.5$, while $\gamma$ is varied [39]. Thus, the log-normal
distribution becomes wider and more skewed as $\gamma$ increases. The number
of fermions is conserved and we considered the $N=L/2$ sector for even $L$ and
$N=(L+1)/2$ sector for odd $L$. Resulting in a Hilbert space size of
$K=\binom{L}{N}$, and a matrix size of $K\times K$.

##  III Short Energy Scales

As a first step we would like to probe the nearest neighbor level spacing
statistics of the CSYK box distributed in order to establish the extended
regime of this model. One expects that in the extended regime short energy
scales (long times) follow the Wigner statistics. The short energy scale
statistics is revealed by the ratio statistics, defined as:

$\displaystyle r$ $\displaystyle=$
$\displaystyle\langle\min(r_{n},r_{n}^{-1})\rangle,$ (3) $\displaystyle r_{n}$
$\displaystyle=$ $\displaystyle\frac{E_{n}-E_{n-1}}{E_{n+1}-E_{n}},$

where $E_{n}$ is the n-th eigenvalue of the Hamiltonian and
$\langle\ldots\rangle$ is an average over an ensemble of different
realizations of disorder and half of the eigenvalues around the middle of the
band. For the Wigner distribution $r_{s}\cong 0.5307$ for the GOE symmetry,
$r_{s}\cong 0.5996$ for GUE and $r_{s}\cong 0.6744$ for GSE. [40].

An interesting behavior emerges for the CSYK. It is known that as a
consequence of the symmetries of the SYK model, the spectrum of Eq. (1) shows
statistics which depends on $L$ [6, 7, 8]. For $L\mod 4=0$ the statistics are
GOE, For $L\mod 4=2$ the statistics are GSE, and for $L\mod 4=1,3$ the
statistics are GUE. This is indeed seen in Fig. 1 where $r_{s}$ averaged over
the middle half of the energy spectrum, for different sizes
$L=8,9,10,11,12,13,14,15,16,17$ and number of fermions $N=4,5,5,6,6,7,7,8,8,9$
is plotted. In all cases we exactly diagonalize the corresponding $K\times K$
matrix and average over $3000$ realizations of disorder (except for the
largest size for which only $100$ realization were computed). It can be seen
that the expected 4-fold symmetry is followed.

Figure 1:  The nearest neighbor level spacing statistics as manifested in the
behavior of the ratio statistics $r_{s}$ for different system sizes $L$ of the
CSYK with a box distribution are indicated by the black circles. The number of
fermions is $N=L/2$ for even $L$ and $N=(L+1)/2$ for odd $L$. The $r_{s}$
values expected for GOE (dashed red) GUE (dashed green) and GSE (dashed blue)
are marked. The expected 4-fold symmetry is seen.

Concentrating on the $L=16$ with $N=8$ systems, we investigate the role played
by changing the distribution from the box distribution to a log-normal.
Setting $p=0.5$ and increasing $\gamma$ we sweep through the values
$\gamma=1,1.5,2,2.5,3,4,5$. For all these values $r_{s}\cong 0.5307\pm 0.001$.
Thus, GOE universal behavior on short energy scales is perfectly followed.

##  IV Local ensemble unfolding

As discussed, the practice of determining the Thouless energy is fraught with
difficulties. In order to compare any spectrum to RMT predictions, one must
recast the spectrum such that it will exhibit an averaged constant density of
states, i.e., a constant level spacing throughout the region examined. What is
the averaging procedure? Usually, one averages the level spacing over an
ensemble of disordered realizations in a given region and then reconstructs a
particular spectrum such as the level spacing is on the average equal to 1
everywhere. Specifically, the averaged level spacing for the i-th level is
$\Delta_{i}=\langle E_{i+m}-E_{i-m}\rangle/2m$ (where $m$ is $O(1)$, here
chosen as $m=5$), and the unfolded spectrum for the j-th realization is
$\varepsilon^{j}_{i}=\varepsilon^{j}_{i-1}+(E^{j}_{i}-E^{j}_{i-1})/\Delta_{i}$.
For brevity we shall call this unfolding procedure local ensemble unfolding.

Implementation of this local unfolding procedure for the level number variance
of CSYK with a box distribution and for a log-normal distribution
$\gamma=1,1.5$ results in the behavior depicted in the inset of Fig 2, which
is in agreement with the behavior observed in Refs. [6, 30]. Here the number
variance begins by following the GOE predictions and grows quadratically for
larger $\langle n\rangle$. The Thouless energy corresponds to the point where
the number variance deviates from GOE predictions. As an estimate of the
Thouless energy we chose the point for which the variance deviates by an
arbitrary amount (set as $0.2$), resulting in $E_{Th}\sim 95\Delta$, for the
box distribution, $E_{Th}\sim 13\Delta$, for $\gamma=1$, and $E_{Th}\sim
7\Delta$ for $\gamma=1.5$, where $\Delta=\langle\Delta_{i}\rangle$. This will
naturally lead to the conclusion that as $\gamma$ increases, $E_{Th}$ strongly
decreases. For larger values of $\gamma$ the variance departs from the GOE
predictions close to $\Delta$ and were not plotted to avoid cluttering the
figure at small $\langle n\rangle$.

Figure 2:  The level number variance, $\langle\delta^{2}n(E)\rangle$, as
function of the energy window $E$. Symbols (black - box, red - $\gamma=1$,
blue $\gamma=1.5$ represent the variance with local ensemble unfolding, fitted
for larger values by $a+b\langle n(E)\rangle^{2}$ (where $b=1.78\ 10^{-5}$,
for the box distribution, $b=1.25\ 10^{-3}$, for $\gamma=1$, and $b=4.4\
10^{-3}$, for $\gamma=1.5$. The cyan line is the GOE prediction
$\langle\delta^{2}n(E)\rangle=(2/\pi^{2})\ln(\langle n(E)\rangle)+0.44$.
Insert: zoom into smaller values of $<n>$. Deviation from the GOE behavior are
observed. A curve depicting GOE plus a constant of $0.2$ corresponds the
dashed cyan line. Using the intersection between the variance and the dashed
cyan line to determine the Thouless energy results in $E_{Th}\sim 95\Delta$
for the box distribution, $E_{Th}\sim 13\Delta$, for $\gamma=1$, and
$E_{Th}\sim 7\Delta$ for $\gamma=1.5$.

As emphasized by Jia and Verbaarschot [30], since unlike typical RMT models
for which all non-diagonal terms are are random, SYK models have just
$\binom{2L}{4}\ll\binom{L}{N}$ independent non-diagonal terms, which leads to
significant sample to sample fluctuations within the ensemble and using the
average level spacing obtained by ensemble average may significantly skew the
number variance. Thus, we need a way to characterize the level spacing for a
specific sample. As previously discussed, in Ref. 30 this was achieved by
parameterized the spectrum using Q-Hermite orthogonal polynomials. We will use
the SVD method as described in the next section.

##  V Singular Value Decomposition

SVD can be used to characterize the features of the spectrum, for example on
what scale does it follow RMT predictions [32, 33, 34, 35, 36, 37, 38]. For
the analysis, one tabulates $P$ eigenvalues around the center of the band of
$M$ realizations of disorder as a matrix $X$ of size $M\times P$ where
$X_{mp}$ is the $p$ level of the $m$-th realization. The matrix $X$ is
decomposed to $X=U\Sigma V^{T}$, where $U$ and $V$ are $M\times M$ and
$P\times P$ matrices correspondingly, and $\Sigma$ is a diagonal matrix of
size $M\times P$ and rank $r=\min(M,P)$. The $r$ diagonal elements of $\Sigma$
are the singular values amplitudes (SV) $\sigma_{k}$ of $X$. All $\sigma_{k}$
are positive and therefore may be ordered by their size
$\sigma_{1}\geq\sigma_{2}\geq\ldots\sigma_{r}$. The Hilbert-Schmidt norm of
the matrix $||X||_{HS}=\sqrt{TrX^{{\dagger}}X}=\sum_{k}\lambda_{k}$ (where
$\lambda_{k}=\sigma_{k}^{2}$). The matrix $X$ could be written as a series
composed of matrices $X^{(k)}$, where $X^{(k)}_{ij}=U_{ik}V^{T}_{jk}$ and
$X_{ij}=\sum_{k}\sigma_{k}X^{(k)}_{ij}$. This series in an approximation of
matrix $X$, where the sum of the first $m$ modes gives a matrix
$\tilde{X}=\sum_{k=1}^{m}\sigma_{k}X^{(k)}$, for which
$||X||_{HS}-||\tilde{X}||_{HS}$ is minimal.

The idea is that low modes capture the long-wavelength fluctuations of the
spectrum, while higher modes sample the shorter wavelength fluctuations. If a
distinct pattern of behavior of the amplitude as function of the mode number
can be seen for a particular range of $k$, it is meaningful to discuss
different behaviors at different energy scales. Indeed, examining the scree
plot of the singular values $\lambda_{k}$ for box distribution and log-normal
distribution with different values of $\gamma$ and a fixed $p=0.5$, for
$L=16$, $N=8$, with $M=4096$ realizations and $P=4096$ eigenvalues around the
middle of the band one sees two distinct regions. The lowest modes ($k=1,2$
for box and $\gamma=1$; $k=1,2,3,4$ for $1.5\leq\gamma\leq 2.5$; and
$k=1,2,3,4,5,6$ for $3\leq\gamma\leq 5$) amplitudes are much larger than all
the other modes, and determine the very long-wave behavior of the spectrum.
The bulk of the modes for larger $k$ follow a power-law behavior
($\lambda_{k}\sim k^{-\alpha}$) with $\alpha=1$, as expected for Wigner-Dyson
statistics [32, 33, 34].

Figure 3:  The scree plot of the singular values for the CSYK with box and
log-normal distribution with $p=0.5$ and $\gamma=1,1.5,2,2.5,3,4,5$ for
$L=16$, $N=8$, with $M=4096$ realizations and $P=4096$ eigenvalues around the
middle of the band. The square amplitude of the singular vale $\lambda_{k}$
are indicated by the symbols. The cyan line corresponds to $\lambda_{k}\sim
k^{-\alpha}$, with $\alpha=1$, as expected for a spectrum which follows
Wigner-Dyson statistics. The lower $k$ modes which capture the non-universal
global structure of the spectrum deviate from Wigner-Dyson.

In order to illustrate the difference between the local ensemble unfolding and
unfolding using the lower modes of the SVD it is useful to examine the
difference in the behavior of level spacing of the i-th level, $\Delta_{i}$
obtained by each method. While for the local ensemble unfolding the level
spacing $\Delta_{i}$ is averaged over all realizations and therefore is not
realization dependent, for the SVD unfolding the i-th level spacing of the
j-th realization
$\Delta^{j}_{i}=(\tilde{\varepsilon}^{j}_{i+m}-\tilde{\varepsilon}^{j}_{i-m})/2m$,
where $\tilde{\varepsilon}^{j}_{i}=\sum_{k=1}^{2}\sigma_{k}U_{ik}V^{T}_{jk}$,
is realization dependent. As can be seen in Fig. 4 there is a noticeable
difference between the realization specific level spacing $\Delta^{j}_{i}$ and
the ensemble average $\Delta_{i}$. This difference becomes larger as $\gamma$
increases, i.e.. as the width of the log-normal distribution increases.
Moreover, the difference, $\Delta^{j}_{i}-\Delta_{i}$, is long-range
correlated for a given realization across thousands of levels. Thus, it makes
sense to define the typical spacing difference between realizations
$\delta\Delta=\sqrt{\langle(\Delta^{j}_{i}-\Delta_{i})^{2}\rangle_{i}}$. For
$L=16$, $N=8$ with $M=4096$ realizations and $P=4096$ levels around the center
of the band one gets $\delta\Delta=4.44\ 10^{-3}\Delta$ for the box
distribution $\delta\Delta=3.52\ 10^{-2}\Delta$ for $\gamma=1$ and
$\delta\Delta=6.61\ 10^{-2}\Delta$ for $\gamma=1.5$.

Figure 4:  The level spacing for the CSYK with box distribution and different
log-normal distribution with $\gamma=1,2,3,4,5$ for $L=16$, $N=8$, for
$M=4096$ realizations and $P=4096$ eigenvalues around the middle of the band.
The ensemble averaged level spacing, $\Delta_{i}$, for the local unfolding
method is the same for all samples and is depicted by the black curve. The SVD
unfolded level spacing , $\Delta^{j}_{i}$, is realization dependent. Five
individual realizations for each distribution are shown
(red,green,blue,yellow,brown curves). As $\gamma$ increases the sample to
sample fluctuations increases. It is clear that the long range correlations
within a sample are very significant.

This long-range sample to sample level fluctuations can explain the behavior
of the level number variance of the local ensemble unfolding seen in Fig. 3.
Essentially, for larger values of $\langle n\rangle$,
$\langle\delta^{2}n\rangle=a+b\langle n\rangle^{2}$ with $b=1.78\ 10^{-5}$,
for the box distribution, $b=1.25\ 10^{-3}$, for $\gamma=1$, and $b=4.4\
10^{-3}$, for $\gamma=1.5$. The quadratic behavior could be understood as the
consequence of realization specific long-range fluctuation of the level
spacing. Calculating $\langle\delta^{2}n\rangle=\langle(n-\langle
n\rangle)^{2}\rangle$, taking into account that after unfolding $\langle
n\rangle=n$ and for a typical realization $n\sim n(1+\delta\Delta/\Delta)$.
Thus, $\langle\delta^{2}n\rangle\sim(n\delta\Delta/\Delta)^{2}$, and after
plugging in the above mentioned values of $\delta\Delta$ one obtains
$\langle\delta^{2}n\rangle\sim 1.97\ 10^{-5}n^{2}$ for the box distribution,
$\langle\delta^{2}n\rangle\sim 1.24\ 10^{-3}n^{2}$ for $\gamma=1$ and
$\langle\delta^{2}n\rangle\sim 4.36\ 10^{-5}n^{2}$ for $\gamma=1.5$ in good
agreement with the numerical value quoted above. Moreover, using our previous
definition of the Thouless energy as the energy for which the deviation from
RMT results becomes larger than some threshold, $b(E_{TH}/\Delta)^{2}=0.2$,
and using $b=(\delta\Delta/\Delta)^{2}$, one gets
$E_{Th}=0.44\Delta^{2}/\delta\Delta$ resulting in
$E_{Th}=100\Delta,18.5\Delta,6.7\Delta$ for box and log-normal distributions
with $\gamma=1,1.5$ correspondingly, in good agreement with the results in
Fig. 3.

One concludes that the main contribution to the local ensemble average
Thouless energy, $E_{Th}$, comes from the sample to sample long range
fluctuations which can be quantified by $\delta\Delta$, the typical level
spacing difference between the different realizations in the ensemble.

Figure 5:  The SVD unfolded level number variance,
$\langle\delta^{2}n(E)\rangle$, as function of the energy window $E$ for box
and log-normal $\gamma=1,1.5,2,2.5,3,4,5$ represented by the symbols. The cyan
dashed lines line is the GOE prediction
$\langle\delta^{2}n(E)\rangle=(2/\pi^{2})\ln(\langle n(E)\rangle)+0.44$. The
SVD unfolded number variance departs from GOE predictions at $E_{Th^{*}}\sim
800\Delta$ for the box and $\gamma=1$ distributions, at $E_{Th^{*}}\sim
500\Delta$ for $1.5\leq\gamma\leq 2.5$ distributions, and at $E_{Th^{*}}\sim
300\Delta$ for $3\leq\gamma\leq 5$ distributions.

Thus, using the sample specific level spacing from the SVD (i.e.,
$\Delta^{j}_{i}$) for unfolding will eliminate the sample to sample
fluctuations contribution to the number variance. Indeed, using the SVD
unfolded spectrum of the $j$-th realization $\tilde{\varepsilon}^{j}_{i}$
defined by
$\varepsilon^{j}_{i}=\varepsilon^{j}_{i-1}+(\tilde{\varepsilon}^{j}_{i}-\tilde{\varepsilon}^{j}_{i-1})/\Delta^{j}_{i}$,
for the calculation of the realization specific level number variance one
obtains the values depicted in Fig. 5. The number variance fits well with GOE
expectations up to $E_{Th^{*}}\sim 800\Delta$ for the box and $\gamma=1$
distributions, $E_{Th^{*}}\sim 500\Delta$ for $1.5\leq\gamma\leq 2.5$
distributions, and $E_{Th^{*}}\sim 300\Delta$ for $3\leq\gamma\leq 5$
distributions. The deviation is towards lower variance than predicted by GOE,
similar to the behavior of the number variance after unfolding using Q-Hermite
orthogonal polynomials [30]. As shown in Ref. 37 the energy of deviation from
the universal behavior can be read off the scree plot by defining the mode for
which the behavior deviates from the $k^{-1}$ power law as $k_{Th^{*}}$. Thus,
as can be garnered from Fig. 3, $k_{Th^{*}}\sim 2$ for the box and $\gamma=1$
distributions, $k_{Th^{*}}\sim 4$ for $1.5\leq\gamma\leq 2.5$ distributions,
and $k_{Th^{*}}\sim 6$ for $3\leq\gamma\leq 5$ distributions. Estimating the
Thouless energy by $E_{Th^{*}}=r\Delta/2k_{Th^{*}}$ [37] results in
$E_{Th^{*}}\sim 1000\Delta$, $E_{Th^{*}}\sim 500\Delta$, and $E_{Th^{*}}\sim
300\Delta$, not to far from the estimations obtained from Fig. 5.

Figure 6:  The level number variance, $\langle\delta^{2}n(E)\rangle$, as
function of the energy window $E$ for box distribution using local ensemble
unfolding (black symbols) and SVD unfolding (red symbols) for $L=15$, $N=8$,
$P=3000$ and $M=3000$. The cyan dashed lines line is the GUE prediction
$\langle\delta^{2}n(E)\rangle=(1/\pi^{2})\ln(\langle n(E)\rangle)+0.35$, while
green dashed line follows $a+bn^{2}$ with $b=2.48\ 10^{-5}$. The different
behavior for the two unfolding methods is clear. Inset: the scree plot of the
singular vale $\lambda_{k}$ indicated by the symbols as function of $k$ . The
purple line corresponds to $\lambda_{k}\sim k^{-\alpha}$, with $\alpha=1$.

Up to now we have mainly presented results for $L=16$ which is the largest
system for which we have ample statistics. Nevertheless, as can be seen in
Fig. 1, for the CSYK the statistics changes for different values of $L$. In
order to verify that our conclusions regarding the non-self-averaging behavior
of the model are relevant for any symmetry, we calculated
$\langle\delta^{2}n(E)\rangle$ for the box distribution using both local
ensemble unfolding and SVD unfolding for $L=15$, with $N=8$ particles, for
$P=3000$ eigenvalues around the middle of the band for $M=3000$ realizations.
For this case GUE behavior, i.e.,
$\langle\delta^{2}n(E)\rangle=(1/\pi^{2})\ln(\langle n(E)\rangle)+0.35$.
Indeed, as can be seen in Fig. 6, both unfolding methods fit the GUE
expectation for small $\langle n\rangle$. On larger scales the same difference
that was seen for $L=16$ GOE is seen for $L=15$ GUE systems. The local
ensemble unfolding results in $E_{Th}\sim 90\Delta$, and a fit to a quadratic
behavior of the form $a+bn^{2}$ leads to $b=2.48\ 10^{-5}$, while for the SVD
unfolding $E_{Th^{*}}\sim 800\Delta$. Comparing with the typical spacing
difference between realizations for the box distribution which for $L=15$
equals to $\delta\Delta=5.17\ 10^{-3}$, one would estimates
$E_{Th}=0.44\Delta^{2}/\delta\Delta\sim 85\Delta$ and
$b=(\delta\Delta/\Delta)^{2}\sim 2.67\ 10^{-5}$, both in good agreement with
the numerical results. As explained above for
$E_{Th^{*}}=r\Delta/2k_{Th^{*}}$, and here $r=3000$ and $k_{Th^{*}}\sim 2$,
resulting in $E_{Th^{*}}\sim 750\Delta$, again in good agreement with the
results.

##  VI Discussion

Much of the fascination with the SYK is connected to its chaotic behavior.
Since short time scales are associated with large energy scales, deviation
from GOE behavior of the spectra on large energy scales indicate non-chaotic
behavior on short time scales. Estimating the time for which the chaotic
behavior emerges is relevant to the estimation of the scrambling time of the
SYK models, which motivates its application to studies of quantum gravity in
black holes [41]. Ensemble averaging also plays an important role in the
duality between SYK and classical general relativity.

For self-averaging systems there is no difference between averaging over the
ensemble or averaging over a large single realization. As discussed, for
finite size SYK samples, there is a huge difference between the Thouless time
for a realization adjusted unfolding compared to the ensemble averaged
unfolding, a difference which is larger as the distribution of the interaction
elements is more skewed.

Nevertheless, it is relevant to extrapolate to infinite systems
($L\rightarrow\infty$) in order to see whether this difference persists. SYK
has three time scales [42]: (i) band structure time $t_{B}$ associated with
the band width, $B$, (ii) Thouless time, $t_{Th}$ or $t_{Th^{*}}$, on which
much of the paper concentrated and (iii) Heisenberg time, $t_{H}$. Since the
band width depends linearly on $L$ [43], $t_{B}=\hbar/B\sim\hbar/L$. Thus, the
Heisenberg time $t_{H}=\hbar/\Delta$ where $\Delta\sim B/2^{L}$, is
$t_{H}=\hbar\ 2^{L}/L$. The Thouless time for the ensemble average is
$t_{Th}=\hbar/E_{Th}$, where $E_{Th}\sim\Delta L^{2}$ and resulting in
$t_{Th}=\hbar\ 2^{L}/L^{3}$. For the realization adjusted unfolding,
$E_{Th^{*}}\sim 2^{L}\Delta B/k_{Th^{*}}$, with $k_{Th^{*}}\sim O(1)$, thus,
$t_{Th^{*}}\sim k_{Th^{*}}\hbar/B\sim k_{Th^{*}}/L$. Hence, $t_{H}>t_{Th}\gg
t_{Th^{*}}\geq t_{B}$, and while the realization adjusted Thouless time
$\lim_{L\rightarrow\infty}t_{Th^{*}}\rightarrow 0$, the ensemble averaged
Thouless time $\lim_{L\rightarrow\infty}t_{Th}\rightarrow\infty$. The
difference between $t_{Th}$ and $t_{Th^{*}}$ increases as the distribution is
more skewed (Figs. 2 and 5).

This leads to the conclusion that one can not determine the behavior of the
SYK model by ensemble average for times shorter than $t_{Th}$ since shorter
times are non-self-averaging. Moreover, these shorter times ($t_{Th^{*}}\sim
t_{B}<t<t_{Th}$) correspond to a parametically large span of times. Therefore,
when one wishes to study the transition from chaotic to localized behavior in
modified SYK models [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] for which
the number of independent random variables remain low, sample to sample
fluctuations are expected to remain important and non-self-averaging at short
times should be considered. In principal, although the transition occurs at
long times nevertheless these effects could obscure the transition point and
influence the perceived nature of the metallic regime. This will be considered
in future studies.

## References

  * [1] S. Sachdev and J. Ye, Phys. Rev. Lett. 70, 3339 (1993).
  * [2] A. Kitaev, Talks at the KITP on April 7th and May 27th (2015).
  * [3] S. Sachdev, Phys. Rev. X 5, 041025 (2015).
  * [4] J. Maldacena and D. Stanford, Phys. Rev. D 94, 106002 (2016).
  * [5] J. Maldacena, International journal of theoretical physics 38, 1113 (1999).
  * [6] A. M. García-García and J. J. M. Verbaarschot, Phys. Rev. D 94, 126010 (2016).
  * [7] Y-Z. You, A. W. W. Ludwig, and C. Xu, Phys. Rev. B 95, 115150 (2017).
  * [8] T. Li, J. Liu, Yuan Xin and Y. Zhou, J. High Energ. Phys. 06, 111 (2017).
  * [9] A. M. García-García, Y. Jia, and J. J. M. Verbaarschot, Phys. Rev. D 97, 106003 (2018).
  * [10] A. M. García-García , B. Loureiro, A. Romero-Bermúdez, and M. Tezuka, Phys. Rev. Lett. 120, 241603 (2018).
  * [11] T. Nosaka, D. Rosa, and J. Yoon, The Thouless time, J. High Energ. Phys. 06, 41 (2018).
  * [12] A. M. García-García and M. Tezuka, Phys. Rev. B 99, 054202 (2019).
  * [13] T. Micklitz, F. Monteiro, A. Altland, Phys. Rev. Lett. 123, 125702 (2019).
  * [14] F. Monteiro, M. Tezuka, A. Altland, D. A. Huse, T. Micklitz, Phys. Rev. Lett. 127, 030601 (2021).
  * [15] F. Monteiro, T. Micklitz, M. Tezuka, and A. Altland, Phys. Rev. Research 3, 013023 (2021).
  * [16] C.-M. Jian, Z. Bi, and C. Xu, Phys. Rev. B 96, 115122 (2017).
  * [17] S.-K. Jian and H. Yao, Phys. Rev. Lett. 119, 206602 (2017).
  * [18] X. Chen, R. Fan, Y. Chen, H. Zhai, and P. Zhang, Phys. Rev. Lett. 119, 207603 (2017).
  * [19] A. Altland, D. Bagrets, and A. Kamenev, Phys. Rev. Lett. 123, 106601 (2019).
  * [20] D. K. Nandy,T. C̃adež, B. Dietz, A. Andreanov, and D. Rosa, arXiv:2206.08599.
  * [21] M. L. Mehta, Random matrices (Acad. Press, New York, 1991), 2nd ed.
  * [22] B. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides and H. B. Shore, Phys. Rev. B. 47, 11487 (1993).
  * [23] T. Guhr, A. Muller-Groeling, H. A. Weidenmuller, Phys. Rep. 299, 190 (1998).
  * [24] Y. Alhassid, Rev. Mod. Phys. 72, 895 (2000).
  * [25] A.D. Mirlin, Phys. Rep. 326, 259 (2000).
  * [26] R. Evers and A.D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * [27] B. Altshuler and B. Shklovskii, Sov. Phys. JETP [Zh. Eksp. Teor. Fiz. 91,220] 64, 127 (1986).
  * [28] J. Sonner and M. Vielma, J. High Energ. Phys. 11 149 (2017).
  * [29] A. Altland and D. Bagrets, Nucl. Phys. B 930, 45 (2018).
  * [30] Y. Jia and J. J.M. Verbaarschot, J. High Energ. Phys. 07 193 (2020).
  * [31] A. Aharony and A. Harris, Phys. Rev. Lett. 77, 3700 (1996).
  * [32] R. Fossion, G. Torres-Vargas and J. C. López-Vieyra, Phys. Rev. E, 88, 060902(R) (2013).
  * [33] G. Torres-Vargas, R. Fossion, C. Tapia-Ignacio and J. C. López-Vieyra, Phys. Rev. E, 96, 012110 (2017).
  * [34] G. Torres-Vargas, J. A. Méndez-Berm´udez, J. C. LópezVieyra and R. Fossion, Phys. Rev. E, 98, 022110 (2018).
  * [35] R. Berkovits, Phys. Rev. B 102, 165140 (2020).
  * [36] R. Berkovits, Phys. Rev. B 104, 054207 (2021).
  * [37] R. Berkovits, Phys. Rev. B 105, 104203 (2022).
  * [38] W.-J. Rao, Phys. Rev. B 105, 054207 (2022).
  * [39] I. M. Khaymovich, V. E. Kravtsov, B. L. Altshuler, and L. B. Ioffe, Phys. Rev. Res. 2, 043346 (2020).
  * [40] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
  * [41] S. Sachdev, arXiv:2205.02285.
  * [42] A. Altland, D. Bagrets, P. Nayak, J. Sonner, and M. Vielma, Phys. Rev. Res. 3, 033259 (2021).
  * [43] J. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. Shenker, D. Stanford, A. Streicher, and M. Tezuka, J. High Energ. Phys. 05, 118 (2017).


# Measuring The Soft Excess Region Size Relative to the Corona in AGN With
NICER

A. Zoghbi Department of Astronomy, University of Maryland, College Park, MD
20742  HEASARC, Code 6601, NASA/GSFC, Greenbelt, MD 20771  CRESST II, NASA
Goddard Space Flight Center, Greenbelt, MD 20771  J. M. Miller Department of
Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

(Accepted for publication in ApJ)

###### Abstract

The soft excess is a significant emission component in the Soft ($<1$ keV)
X-ray spectra of many AGN. It has been explained by disk reflection, a warm
corona and other models. Understanding its origin is crucial for the energy
budget of AGN emission, and for using it to study the inner accretion disk.
Here, we track the weeks-to-months variability of several AGN that show
different levels of soft excess strength with NICER. We use the variability
time scales to compare the relative size of the soft excess emission region to
the corona producing the hard X-ray emission above 1 keV. We find that the
size of the soft excess emission region relative to the corona is not the same
for the three sources studied. For TON S180, the soft excess region is
comparable in size to the hard corona. While for MRK 335 and 1H0707-495, the
soft excess region is larger than the corona by a factor of 2-4. This is the
first time the relative sizes are quantified independently of the assumptions
of the spectral models.

X-ray active galactic nuclei (2035), Active galactic nuclei (16), Seyfert
galaxies (1447), Supermassive black holes (1663), Black hole physics (159)

††facilities: NICER(XTI), HEASARC

##  1 Introduction

A puzzling feature the X-ray spectra of AGN is the so-called soft excess. This
is a strong excess of emission that is observed above the extrapolation of the
power-law from the hot corona to energies below 1 keV (e.g. Walter & Fink,
1993; Bianchi et al., 2009). The origin of this feature is still debated
(García et al., 2019; Petrucci et al., 2020). It is featureless, with a shape
that can generally be described by a blackbody with a temperature of $\sim$0.5
keV (Gierliński & Done, 2004).

Early models attributed the soft excess to the tail of the disk emission
Leighly (1999) or to a smeared blend absorption lines (Gierliński & Done,
2004). In recent years, two models are often discussed in the literature. In
the first, the excess is produced by the sum of recombination lines and
bremsstrahlung from the heated surface of the disk that is illuminated by the
hot corona. When produced close enough to the black hole, relativistic effects
smear and broaden the emission lines so it appears featureless (Crummy et al.,
2006; Walton et al., 2013; García et al., 2019). Reverberation time lags
observed in many sources are a natural consequence of this model (Zoghbi et
al., 2010; De Marco et al., 2013).

In the second model, the soft excess is produced by Comptonization of thermal
disk photons by a warm ($\sim 1$ keV) and optically thick ($\tau\sim 10-20$)
layer of gas at the surface of the disk (Magdziarz et al., 1998; Czerny et
al., 2003; Done et al., 2012; Petrucci et al., 2013). The warm corona is
distinct from the hot corona producing the primary hard X-ray emission, but
they must be heated by accretion power. Both these models describe the
observed time average spectra equally well, with the main discussion focusing
on whether they are physically consistent (García et al., 2019; Petrucci et
al., 2020; Ballantyne & Xiang, 2020).

These models predict different sizes for the soft excess region. In this work,
we report on a monitoring experiment with NICER that tracks the variability of
the soft excess in several sources, and use the variability to infer the size
of the soft excess region relative to the corona that emits at hard energies
($>2$ keV). The combination of the large effective area and monitoring
capability of NICER allow for this experiment to be conducted for the first
time.

##  2 Analysis

###  2.1 Data

The long term variability of 4 Narrow Line Seyfert 1 galaxies is presented.
Two of which were part of the original NICER proposal (TON S180 and PG
1404+226). Two others (1H0707-495 and MRK 335) that had public data on the
NICER archive were also included.

Our analysis includes all public NICER observations available on Dec 2022. We
use analysis tools available in heasoft v6.31.1. Cleaned events files were
generated from the unfiltered data using nicerl2. The spectra and the
corresponding response and area files were then generated using nicerl3-spect.
Background spectra were generated for every observation using all three
background estimators. We find that the Scorpeon and the 3C50 models were
consistent, while the space weather model over-estimates the background. We
report the results from using the Scorpeon model. Using the 3C50 model gives
similar resul.

Figure 1: Combined NICER spectra from all observations of each sources. These
are $EF(E)$ spectra that have been unfolded against a constant model to factor
out the effective area of the detector. Note that this unfolding leaves some
instrumental features around 2 keV.

The net exposures of the resulting spectra are typically in the range of 0.6–3
ks. 1H0707-495 had a handful of observations with exposures $>10$ ks. The
number of observations for each source are: 20 (1H0707-495), 80 (MRK 335), 29
(TON S180), and 26 (PG 1404+226). To visualize the spectral shape, Figure 1
shows the combined spectra for each of the 4 objects, produced by combining
the individual spectra using the addspec tool. These are $EF(E)$ spectra that
have been unfolded against a constant model to factor out the effective area
of the detector. The fact that the instrument response is averaged over many
months when combining the spectra results in some spurious instrument features
around 2 keV. These do not affect our analysis because we are interested in
the flux from broad spectral components. Figure 2 shows the total (0.3–10 keV)
count rate light curves.

Figure 2: The count rate light curves in the 0.3–10 keV band from all
observations. The x-axes is in units of modified Julian date, and they all
span the same $\sim 580$ days.

Figure 1 shows that the 4 objects span a range in both flux and strength of
the soft excess component. Figure 2 shows that our targeted program for TON
S180 and PG 1404+226 span about a year and with some uniform sampling compared
to 1H0707-495. MRK 335 has the largest number of observations.

In order to characterize the long term variability of the soft excess relative
to the hard component from the corona dominating above 2 keV, we model the
spectra from the individual observations with a model consisting of a power-
law and a disk blackbody component. The power-law models the spectrum of the
corona, dominating above 2 keV, while the disk blackbody models the soft
excess. This model is used because first, we are only interested in measuring
the fluxes of the two components, and second, the 0.8–3 keV does not allow for
more complex models. We find that this modeling is sufficient for the purpose
of measuring the fluxes needed in this work. Using other models for the soft
excess have little effect on the final results.

We note that although PG 1404+226 shows a strong soft excess relative the hard
component (Figure 1), the latter is too faint to enable spectral modeling of
the individual spectra above 2 keV. We therefore drop it from subsequent
analysis that compares the variability of the soft excess component to the
hard corona.

For some observations of 1H0707-495 and MRK 335, the signal above 2 keV was
only high enough to constrain the flux of the power-law component, but not the
spectral index. So for these two sources, we fix the spectral index of the
power-law to the value measured from fitting the total spectrum with a similar
model. For TON S180, the signal is high enough to allow for both the flux and
spectral index to be measured in the individual spectra.

After the spectral modeling of each observation, we obtain a light curve for
the 0.3–10 keV fluxes for the blackbody ($F_{bb}$) and power-law $F_{p}$
components, the blackbody temperature, and for the case of TON S180, for the
photon index too. We then proceed by calculating and modeling the power
spectral density (PSD) from these light curves.

###  2.2 Power Spectra

The light curves are not evenly-sampled. We use the likelihood estimation
method from Zoghbi et al. (2013). We specifically use the fqlag package v0.3.4
(Zoghbi, 2023a). We first estimate the power values at 10 logarithmically-
spaced frequency bins in the observed frequency range. The lowest and highest
frequency limit are $1/T$ and $0.5/\Delta T_{\rm min}$, where $T$ is the
length of the observing campaign, and $\Delta T_{\rm min}$ is the minimum time
separation between neighboring observations. Also, as recommended by package
documentation, we include two buffer frequency bins at both end of these
limits to minimize biases, which are ignored when reporting the PSD.

Figure 3: Results of the PSD modeling for the fluxes of the blackbody
($F_{bb}$) and the coronal power-law components ($F_{p}$) for 1H0707-495, MRK
335 and TON S180. _T_ op: The measured PSD values at 10 frequency bins as
calculated using fqlag. The values and errors are the median and standard
deviation from the MCMC chains. The y-axis is in units of RMS-normalized PSD
and x-axis is the logarithm of the the Fourier frequency. _B_ ottom: Smoothed
density estimates for the probability density of the key parameter of
interest, resulting from modeling the PSDs in the top panel by either a a
zero-centered Lorentzian model (1H0707-495) or power-law model (MRK 335 and
TON S180). For the former case, the logarithm of the break frequency is
plotted, while for the latter, the integrated variability power (integral
under the power-law PSD) in percentage units is plotted.

Because the frequency bins may not be large enough to ensure the errors in the
measured PSD values are always Gaussian, and in order to allow for subsequent
fitting of the PSD, we do not just use a single value estimate of the PSD
(e.g. median and a standard deviation), we instead run Monte Carlo Markov
Chains (MCMC) and empirically measure the probability density of the PSD
values at every frequency bin. These probability density estimates are then
approximated by a flexible standard probability density function (PDF). After
trying several general functions, we found that the Johnson SB distribution
provides an excellent approximation to the MCMC distributions. All the values
were individually inspected to ensure the approximation is adequate.

To characterize the measured PSDs, we fit them with three models: a power-law,
a bending power-law and zero centered Lorentzian. The bending power-law has an
index that smoothly changes between two values at some break frequency
(McHardy et al., 2004). We tested fixing the lower index at both 0 and 1.
Neither of them provided a significant improvement over the other two models
(power-law and zero-centered Lorentzian).

In the end, we found that the power-law model was sufficient to describe the
PSD of both $F_{p}$ and $F_{bb}$ for both MRK 335 and TON S180, while the
zero-centered Lorentzian provided a significantly better fit for 1H0707-495.
In other words, a frequency break is significantly detected only on
1H0707-495. A plot of the PSD shapes in the three cases are shown in the top
panel of Figure 3.

We note that if we use the relation between black hole mass and variability
time scale from McHardy et al. 2006, and mass and Luminosity estimates from
the literature for MRK 335 and TON S180, we estimate the break frequencies to
be: $log\nu\sim$-0.8 and 1.4 $\rm{days}^{-1}$, respectively, which are higher
than our cadence allows. For MRK 335, there is a hint of a break at the top of
Figure 3, consistent with the expected value, but it is not significant.

| 1H0707-495 | MRK 335 | TON S180
---|---|---|---
$F_{p}$ | $A$ | $2.57\pm 0.17$ | $-0.87\pm 0.16$ | $-2.26\pm 0.30$
$\nu_{b}$ | $-0.41\pm 0.19$ | - | -
$\alpha$ | - | $-0.76\pm 0.11$ | $-1.22\pm 0.25$
$F_{bb}$ | $A$ | $1.96\pm 0.13$ | $-1.37\pm 0.15$ | $-2.19\pm 0.29$
$\nu_{b}$ | $-0.99\pm 0.16$ | - | -
$\alpha$ | - | $-0.90\pm 0.15$ | $-1.23\pm 0.20$
Table 1: Summary of the PSD fit parameters. For 1H0707-495, a zero-centered
Lorentzian is fitted to the PSD, so the amplitude $A$ (in units of RMS
$\times$ days) and break frequency $\nu_{b}$ (in units of log
$\rm{days}^{-1}$) are reported. For MRK 335 and TON S180, a power-law model is
fitted, so the amplitude $A$ and index $\alpha$ are reported.

In the bottom of Figure 3, we plot the probability distribution of a summary
parameter that characterizes the variability of the flux of the two
components. In the case of the 1H0707-495, we have a direct measure of the
break frequency, so we use it. For MRK 335 and TON S180, no break in the PSD
is measured, so we use the total integrated RMS variability (i.e. the integral
under the PSD Vaughan et al., 2003). For each source, we compare the
variability of the blackbody ($F_{bb}$) and power-law ($F_{p}$) spectral
components.

Figure 4: Probability density of the relative size of the soft excess to the
hot corona. The relative size is obtained as the ratio of the summary
parameter for soft excess flux ($F_{bb}$) and the corona ($F_{p}$). For
1H0707-495, the summary parameters is the PSD break frequency, while for the
TON S180 and MRK 335, we use the total integrated RMS variability. The
$1\sigma$ single value estimates from these probability densities are:
$3.7_{-1.5}^{+2.7}$, $2.2_{-0.6}^{+1.0}$ and $0.9_{-0.4}^{+1.0}$ for
1H0707-495, MRK 335 and TON S180, respectively.

##  3 Results & Discussion

Section 2.2 (summarized in Figure 3), presented one way to characterize the
variability from the fluxes of the soft excess component, modeled with a
blackbody, and the hard coronal component, modeled with a power-law. It is not
straight forward to convert these measurement to physical units. However, we
can focus on the _relative_ scale of the parameters for the two spectral
components.

Figure 4 shows the ratio between the summary parameters for the soft excess
component ($F_{bb}$) and the corona component ($F_{p}$). The summary parameter
is the characteristic time scale (i.e. the break frequency) for the case
1H0707-495, and the total integrated variability power in the case of the MRK
335 and TON S180. The $1\sigma$ single value estimates from these probability
densities are: $3.7_{-1.5}^{+2.7}$, $2.2_{-0.6}^{+1.0}$ and
$0.9_{-0.4}^{+1.0}$ for 1H0707-495, MRK 335 and TON S180, respectively.

For the case of 1H0707-495, both the break frequency and the total variability
power can be measured. To check for consistency, we also measure the relative
size using the variability power, we find: $4.1^{+2.3}_{-1.5}$, which is
similar to the value measured from the break frequency, supporting the
robustness of the result.

With the assumption that the parameter characterizing the variability (either
the break frequency or the integrated RMS power) scale with the size of the
emission region, the ratio in Figure 4 maps directly to the relative size of
the emission regions producing the soft excess and the hard corona. This
assumption is reasonable, and is justified by the many scaling relations
observed in accreting black holes. This includes the scaling of variability
time-scale with black hole mass (McHardy et al., 2006), the normalized excess
variance ($\sigma_{\rm NXV}^{2}$) scaling with black hole mass (Papadakis,
2004; Ponti et al., 2012; Akylas et al., 2022), and the RMS variability on
long time scales in optical light curves also scaling with mass (MacLeod et
al., 2010).

Note that what we refer to as size here is the radial location of the emission
region, which is equivalent to region scale size given the usual assumption of
symmetry in accretion disks.

The plot in Figure 4 shows that the soft excess region is comparable in size
to the hard corona in TON S180, while it is larger by a factor of 2–4 for the
case of 1H0707-495 and MRK 335. This estimate does not make any assumptions
about the nature of the soft excess emission, and so they provide new insight
into the nature of the soft excess.

Figure 5: Illustration of the different configurations that correspond to the
different models proposed for the soft excess. The left scenario is for a
relativistic reflection model, where the hot corona (blue) illuminates the
disk producing emission lines, which when produced very close to the black
hole are broadened and smeared out. The other two images are for the case of
warm corona: A layer above the disk that is heated by internal dissipation
that is distinct from the hot corona (Ponti et al., 2013), and a co-located
hot and warm corona that radiate inward of some radius $R_{\rm corona}$, while
a standard disk emits outside it (Done et al., 2012).

A key requirement for the relativistic reflection model to explain the smooth
soft excess, is that the emission has to originate very close to the black
hole (only a few gravitation radii; $r_{g}$ from the horizon). This is needed
to blur out emission lines from O and Ne that are produced below 1 keV in in a
partially-ionized gas. This in turn requires the illuminating corona to be
very compact (Fabian et al., 2009; Zoghbi et al., 2010; Wilkins & Fabian,
2012; Walton et al., 2013; Jiang et al., 2019). So the primary X-ray source
has to be very compact (only a few $r_{g}$), and the reflection is also very
compact (Illustrated in the left panel of Figure 5). So to a first order, the
reflector is expected to be _comparable_ in size to the primary corona. We
confirmed this by running ray-tracing simulations around a spinning black
hole. Although, we leave the detailed modeling to future work, our initial
simulations show that for a standard lamp-post, the relative size is undefined
because the corona is assumed to be a point source. If the an extension is
added to the corona, we find that the size of the reflection region (as
measured for example by the radius that encompasses 90% of the emission) is
comparable to the size of the corona, consistent with emissivity profile
studies (Wilkins & Fabian, 2012).

The geometry for the warm corona model may not be as well defined, and it is
typically an assumption in the model. In some models (e.g. Petrucci et al.,
2013) (Illustrated in the middle panel of Figure 5), an optically-thin hot
corona ($kT\sim 100$ keV; $\tau\sim 1$) is present in the inner parts of the
accretion flow, producing the hard emission. The outer accretion flow is a
vertically structured accretion disk, with cold and optically thick matter in
the deeper layers, while the upper layers are composed of an optically-thick
warm corona ($kT\sim 1$ keV; $\tau\sim 10-20$) that is powered by internal
heating. The prediction from these geometries is that the soft excess region
is _larger_ in size than the hot corona.

Other models (e.g. Done et al., 2012) assume that the emission thermalizes to
a (color temperature corrected) blackbody only at large radii (Illustrated in
the right panel of Figure 5). At smaller radii the gravitational energy is
split between powering optically thick Comptonized disk emission, forming the
soft X-ray excess, and an optically thin corona above the disk, forming the
hot corona tail at higher energies. Similar to the reflection case, the
prediction here is that the soft excess region is comparable in size to the
hot corona.

Other variants of the model (Kubota & Done, 2018) assume an inner hot corona
up to $R_{\rm hot}$, then warm corona up to $R_{\rm warm}$ and then the
standard cold disk. In that specific study, a model for the spectral energy
distribution can be fitted to the observations, and the different radii can be
inferred. Those result suggest a relative size of $R_{\rm warm}/R_{\rm
hot}\sim 2-4$ for the sources in that study, which is similar to what we
measure for 1H0707-495 and MRK 335. The geometry is slightly different from
what we sketched in Figure 5, but it would fall under the middle sketch where
the warm corona is outside the hot corona.

Our results suggest that for TON S180, where the size of the soft excess
region is comparable to the hot corona, the relativistic reflection model or
the warm corona that sandwiches a hot corona are possible geometries. On the
other hand, for 1H0707-495 and MRK 335, a warm corona that is larger than the
hot corona appears to be more consistent with the data.

In our discussion of the relative sizes, we are ignoring the detailed
geometries (e.g. spherical vs flat) and viewing angle. For a simple spherical
emission regions, there is only scale size that controls the variability. For
other shapes, say a torus, there are different scales in different directions,
but the variability should be dominated by the largest scale, regardless of
our viewing angle. Consequently, the flux may depend on the angle, but the
variability timescale does not.

We plan to continue monitoring these and other sources with a soft excess to
obtain further constraints.

The code used to produce these results is available on Zenodo (Zoghbi, 2023b).
The data products are available on the Open Science Framework site (Zoghbi,
2023c).

The material is based upon work supported by NASA under award numbers
80GSFC21M0002 and 80NSSC23K0333.

## References

  * Akylas et al. (2022) Akylas, A., Papadakis, I., & Georgakakis, A. 2022, aap, 666, A127, doi: 10.1051/0004-6361/202244162
  * Ballantyne & Xiang (2020) Ballantyne, D. R., & Xiang, X. 2020, mnras, 496, 4255, doi: 10.1093/mnras/staa1866
  * Bianchi et al. (2009) Bianchi, S., Guainazzi, M., Matt, G., Fonseca Bonilla, N., & Ponti, G. 2009, aap, 495, 421, doi: 10.1051/0004-6361:200810620
  * Crummy et al. (2006) Crummy, J., Fabian, A. C., Gallo, L., & Ross, R. R. 2006, mnras, 365, 1067, doi: 10.1111/j.1365-2966.2005.09844.x
  * Czerny et al. (2003) Czerny, B., Nikołajuk, M., Różańska, A., et al. 2003, aap, 412, 317, doi: 10.1051/0004-6361:20031441
  * De Marco et al. (2013) De Marco, B., Ponti, G., Cappi, M., et al. 2013, mnras, 431, 2441, doi: 10.1093/mnras/stt339
  * Done et al. (2012) Done, C., Davis, S. W., Jin, C., Blaes, O., & Ward, M. 2012, mnras, 420, 1848, doi: 10.1111/j.1365-2966.2011.19779.x
  * Fabian et al. (2009) Fabian, A. C., Zoghbi, A., Ross, R. R., et al. 2009, nat, 459, 540, doi: 10.1038/nature08007
  * García et al. (2019) García, J. A., Kara, E., Walton, D., et al. 2019, apj, 871, 88, doi: 10.3847/1538-4357/aaf739
  * Gierliński & Done (2004) Gierliński, M., & Done, C. 2004, mnras, 349, L7, doi: 10.1111/j.1365-2966.2004.07687.x
  * Jiang et al. (2019) Jiang, J., Fabian, A. C., Dauser, T., et al. 2019, mnras, 489, 3436, doi: 10.1093/mnras/stz2326
  * Kubota & Done (2018) Kubota, A., & Done, C. 2018, mnras, 480, 1247, doi: 10.1093/mnras/sty1890
  * Leighly (1999) Leighly, K. M. 1999, apjs, 125, 317, doi: 10.1086/313287
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, apj, 721, 1014, doi: 10.1088/0004-637X/721/2/1014
  * Magdziarz et al. (1998) Magdziarz, P., Blaes, O. M., Zdziarski, A. A., Johnson, W. N., & Smith, D. A. 1998, mnras, 301, 179, doi: 10.1046/j.1365-8711.1998.02015.x
  * McHardy et al. (2006) McHardy, I. M., Koerding, E., Knigge, C., Uttley, P., & Fender, R. P. 2006, nat, 444, 730, doi: 10.1038/nature05389
  * McHardy et al. (2004) McHardy, I. M., Papadakis, I. E., Uttley, P., Page, M. J., & Mason, K. O. 2004, mnras, 348, 783, doi: 10.1111/j.1365-2966.2004.07376.x
  * Papadakis (2004) Papadakis, I. E. 2004, mnras, 348, 207, doi: 10.1111/j.1365-2966.2004.07351.x
  * Petrucci et al. (2013) Petrucci, P. O., Paltani, S., Malzac, J., et al. 2013, aap, 549, A73, doi: 10.1051/0004-6361/201219956
  * Petrucci et al. (2020) Petrucci, P. O., Gronkiewicz, D., Rozanska, A., et al. 2020, aap, 634, A85, doi: 10.1051/0004-6361/201937011
  * Ponti et al. (2012) Ponti, G., Papadakis, I., Bianchi, S., et al. 2012, aap, 542, A83, doi: 10.1051/0004-6361/201118326
  * Ponti et al. (2013) Ponti, G., Cappi, M., Costantini, E., et al. 2013, aap, 549, A72, doi: 10.1051/0004-6361/201219450
  * Vaughan et al. (2003) Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, mnras, 345, 1271, doi: 10.1046/j.1365-2966.2003.07042.x
  * Walter & Fink (1993) Walter, R., & Fink, H. H. 1993, aap, 274, 105
  * Walton et al. (2013) Walton, D. J., Nardini, E., Fabian, A. C., Gallo, L. C., & Reis, R. C. 2013, mnras, 428, 2901, doi: 10.1093/mnras/sts227
  * Wilkins & Fabian (2012) Wilkins, D. R., & Fabian, A. C. 2012, mnras, 424, 1284, doi: 10.1111/j.1365-2966.2012.21308.x
  * Zoghbi (2023a) Zoghbi, A. 2023a, zoghbi-a/fqlag: v0.3.4, v0.3.4, Zenodo, doi: 10.5281/zenodo.8309628
  * Zoghbi (2023b) —. 2023b, Code for Article: ”Measuring The Soft Excess Region Size Relative to the Corona in AGN With NICER”, 1.0, Zenodo, doi: 10.5281/zenodo.8310123
  * Zoghbi (2023c) —. 2023c, Data Products for “Measuring The Soft Excess Region Size Relative to the Corona in AGN With NICER”, OSF, doi: 10.17605/OSF.IO/N35FT
  * Zoghbi et al. (2010) Zoghbi, A., Fabian, A. C., Uttley, P., et al. 2010, mnras, 401, 2419, doi: 10.1111/j.1365-2966.2009.15816.x
  * Zoghbi et al. (2013) Zoghbi, A., Reynolds, C., & Cackett, E. M. 2013, apj, 777, 24, doi: 10.1088/0004-637X/777/1/24