# Eigenstate thermalization scaling in approaching the classical limit

Goran Nakerst  Department of Theoretical Physics, Maynooth University, Co.
Kildare, Ireland  Masudul Haque  Department of Theoretical Physics, Maynooth
University, Co. Kildare, Ireland  Max-Planck-Institut für Physik komplexer
Systeme, D-01187 Dresden, Germany

###### Abstract

According to the eigenstate thermalization hypothesis (ETH), the eigenstate-
to-eigenstate fluctuations of expectation values of local observables should
decrease with increasing system size. In approaching the thermodynamic limit —
the number of sites and the particle number increasing at the same rate — the
fluctuations should scale as $\sim D^{-1/2}$ with the Hilbert space dimension
$D$. Here, we study a different limit — the classical or semiclassical limit —
by increasing the particle number in fixed lattice topologies. We focus on the
paradigmatic Bose-Hubbard system, which is quantum-chaotic for large lattices
and shows mixed behavior for small lattices. We derive expressions for the
expected scaling, assuming ideal eigenstates having Gaussian-distributed
random components. We show numerically that, for larger lattices, ETH scaling
of physical mid-spectrum eigenstates follows the ideal (Gaussian) expectation,
but for smaller lattices, the scaling occurs via a different exponent. We
examine several plausible mechanisms for this anomalous scaling.

##  I Introduction

During the last decade and half, a considerable amount of research has focused
on understanding how isolated quantum systems can relax and thermalize. A
cornerstone of this understanding is the Eigenstate Thermalization Hypothesis
(ETH) Deutsch (1991); Srednicki (1994, 1996, 1999); Rigol  _et al._ (2008);
Polkovnikov  _et al._ (2011); Rigol and Srednicki (2012); Eisert  _et al._
(2015); D’Alessio  _et al._ (2016). According to the ETH, the diagonal matrix
elements of observables in the eigenstate basis of the Hamiltonian, the
eigenstate expectation values (EEV), coincide locally with the micro-canonical
ensemble. This means that the EEVs vary smoothly as a function of energy
eigenvalues. The smoothness is quantified as the fluctuations of the EEVs
being exponentially small as a function of the system size.

A standard quantitative statement of ETH states that the matrix elements of an
operator $A$ representing a typical physical observable should have the form
Srednicki (1996, 1999)

$\langle
E_{\alpha}|A|E_{\beta}\rangle=\delta_{\alpha\beta}f_{A}^{(1)}(\bar{E})+e^{-S(\bar{E})/2}f_{A}^{(2)}(\bar{E},\omega)R_{\alpha\beta}$
(1)

where $S$ is the entropy, $|E_{\alpha}\rangle$ is an energy eigenstate with
eigenenergy $E_{\alpha}$, $\bar{E}=(E_{\alpha}+E_{\beta})/2$, and
$\omega=E_{\beta}-E_{\alpha}$. The $f_{A}^{(1/2)}$ are smooth functions, and
$R_{\alpha\beta}$ is a (pseudo) random variable with zero mean and unit
variance. For systems with finite Hilbert space dimension $D$, the entropy
scales as $S\sim\log D$. Thus a crucial aspect of ETH is the scaling of the
width of the distribution of either diagonal or off-diagonal matrix elements:
when approaching the thermodynamic limit, this width falls off as
$e^{-S/2}\sim D^{-1/2}$, i.e., exponentially with system size. This scaling
can be understood using the similarity between typical many-body eigenstates
and random states Neuenhahn and Marquardt (2012); Beugeling  _et al._ (2014,
2015a). This behavior contrasts sharply with integrable systems, which do not
obey ETH scaling — the width of diagonal matrix element distributions
generally have power law decay with system size Ziraldo and Santoro (2013);
Ikeda  _et al._ (2013); Beugeling  _et al._ (2014); Alba (2015); Nandy  _et
al._ (2016); Magán (2016); Haque and McClarty (2019); Mierzejewski and Vidmar
(2020), and the off-diagonal matrix element generally has a non-gaussian
distribution Beugeling  _et al._ (2015a); Haque and McClarty (2019); LeBlond
_et al._ (2019).

Evidence from a large number of numerical studies strongly suggests that ETH
is satisfied for typical states of generic nonintegrable systems and for
physical observables Rigol  _et al._ (2008); Rigol (2009); Biroli  _et al._
(2010); Rigol and Santos (2010); Santos and Rigol (2010); Roux (2010);
Motohashi (2011); Neuenhahn and Marquardt (2012); Brandino  _et al._ (2012);
Steinigeweg  _et al._ (2013); Kim  _et al._ (2014); Beugeling  _et al._
(2014); Sorg  _et al._ (2014); Steinigeweg  _et al._ (2014); Beugeling  _et
al._ (2015a); Fratus and Srednicki (2015); Nandkishore and Huse (2015); Khodja
_et al._ (2015); D’Alessio  _et al._ (2016); Mondaini  _et al._ (2016);
Chandran  _et al._ (2016); Luitz and Bar Lev (2016); Mondaini and Rigol
(2017); Sonner and Vielma (2017); Lan and Powell (2017); Yoshizawa  _et al._
(2018); Dymarsky  _et al._ (2018); Hunter-Jones  _et al._ (2018); Haque and
McClarty (2019); Jansen  _et al._ (2019); Khaymovich  _et al._ (2019);
Mierzejewski and Vidmar (2020). In particular, several studies have examined
the decrease of the width of matrix element distributions as the thermodynamic
limit is approached, both for chaotic and integrable systems Steinigeweg  _et
al._ (2013); Ziraldo and Santoro (2013); Ikeda  _et al._ (2013); Beugeling
_et al._ (2014); Steinigeweg  _et al._ (2014); Beugeling  _et al._ (2015a);
Alba (2015); Nandy  _et al._ (2016); Magán (2016); Chandran  _et al._ (2016);
Lan and Powell (2017); Hunter-Jones  _et al._ (2018); Hamazaki and Ueda
(2018); Khemani  _et al._ (2019); Haque and McClarty (2019); Jansen  _et al._
(2019); Khaymovich  _et al._ (2019); Mierzejewski and Vidmar (2020); Brenes
_et al._ (2020). For lattice systems, the thermodynamic limit involves
increasing both the lattice size and the particle number simultaneously,
keeping the average density fixed. In the present work, we will instead
explore approaching the _classical limit_ — we consider increasing particle
numbers in fixed lattice topologies. We investigate the scaling behavior of
EEV fluctuations as a function of Hilbert space dimension, as the classical
limit is approached.

We focus on bosonic systems and study the Bose-Hubbard Hamiltonian on fixed
numbers of sites. The dynamics of these systems in the classical limit have
been extensively studied, where they are described by a discrete nonlinear
Schroedinger equation. Because of the existence of a classical limit, Bose-
Hubbard systems have also been extensively used as a testbed for semi-
classical methods Polkovnikov (2003); Mahmud  _et al._ (2005); Hiller  _et
al._ (2006); Mossmann and Jung (2006); Graefe and Korsch (2007); Trimborn  _et
al._ (2009); Polkovnikov (2010); Pawłowski  _et al._ (2011); Simon and Strunz
(2012, 2014); Veksler and Fishman (2015); Engl  _et al._ (2015, 2016); Ray
_et al._ (2016a); Dubertrand and Müller (2016); Tomsovic  _et al._ (2018);
Tomsovic (2018); Rammensee  _et al._ (2018); Kidd  _et al._ (2019); Schlagheck
_et al._ (2019). Quantum dynamics of these systems have also been compared
with the dynamics of the corresponding classical limit Nemoto  _et al._
(2000); Weiss and Teichmann (2008); Viscondi and Furuya (2011); Gertjerenken
and Weiss (2013); Kolovsky (2016); Heinisch and Holthaus (2016); Rautenberg
and Gärttner (2020).

We consider $N$ particles in $k$ sites, arranged linearly. We approach the
large limit of large Hilbert space dimension ($D\to\infty$) by keeping $k$
fixed and increasing $N$, as opposed to the usual thermodynamic limit for
which the ratio $N/k$ would be kept fixed. The 2-site system (Bose-Hubbard
dimer, $k=2$) is integrable, and hence we omit this case in the present work
on ETH scaling. We focus on lattices with sizes from $k=3$ up to $k=10$. The
case of $k=3$ (and to a lesser extent $k=4$) is particularly interesting:
although not integrable, the classical phase space in this case is known to be
highly ‘mixed’ Mossmann and Jung (2006); Hiller  _et al._ (2009); Viscondi and
Furuya (2011); Arwas  _et al._ (2014); Han and Wu (2016); Bürkle and Anglin
(2019); Bychek  _et al._ (2020); Rautenberg and Gärttner (2020), and the
behavior of the finite-$N$ quantum system shows deviations from fully chaotic
behavior. For larger site numbers, $k\gtrsim 5$, the Bose-Hubbard systems tend
to behave as quantum-chaotic systems Kolovsky and Buchleitner (2004); Kollath
_et al._ (2010); Kolovsky (2016); de la Cruz  _et al._ (2020); Pausch  _et
al._ (2020).

Assuming the eigenstates to be perfectly ergodic “infinite-temperature” states
in the sense of having Gaussian-distributed random coefficients, an exact
formula for the EEV fluctuations is available. This relation [Eq. (7)] only
depends on trace expressions of operators, and holds for all operators, as
long as the “Gaussian eigenstates” approximation is valid. It mathematically
connects the Gaussianity of eigenstates to finite size scaling of ETH, and
provides an alternate derivation of ETH scaling.

The class of observables generally considered in ETH investigations
(observables of physical interest) are few-body operators or sums thereof,
i.e., ‘local’ operators. However, having $k$ fixed and small means that the
notion of ‘local’ operators has to be re-examined. When considering
$k\to\infty$, two-site operators are local in the sense that $2\ll{k}$. This
is no longer true for our constant $k$; in particular it is strongly violated
for $k=3$. Therefore we re-examine the expected scaling of EEV fluctuations,
first for Gaussian eigenstates. We show that the expected behavior for
$N\gg{k}$ is $\sim D^{-e_{0}}$ with $e_{0}=\frac{1}{2}-\frac{1}{k-1}$. The
exponent ranges from $e_{0}=0$ for the 3-site case ($k=3$) to
$e_{0}\to\frac{1}{2}$ for large $k$ ($N\gg{k}\gg 1$).

The scaling exponent $e_{0}$ being different from 1/2 for moderate $k$ can be
traced to the unboundedness of the operators in the $N\to\infty$ limit. The
$D^{-1/2}$ scaling could be retrieved by rescaling the operators $A$ and
considering operators $\bar{A}=A/N$.

For the actual mid-spectrum eigenstates of the Bose-Hubbard chains, we show
numerically that the EEV fluctuations match the behavior predicted by the
Gaussian ansatz, once $k$ is large enough. For small $k$, the dependence on
$D$ is seen to be a power law but with a _different_ exponent $e\neq{e_{0}}$.
This is a remarkable manifestation of the ‘mixed’ nature of the 3-site and
4-site Bose-Hubbard systems.

At present, we do not have a predictive explanation for the observed
exponents. Small or subleading deviations from Gaussian behavior (from random-
matrix behavior) has been noted previously in mid-spectrum eigenstates of
chaotic many-body systems. Such deviations can be seen in the eigenstate
coefficient distribution in the configuration basis Beugeling  _et al._
(2018); Luitz and Bar Lev (2016); Bäcker  _et al._ (2019); Luitz  _et al._
(2020), in the entanglement entropy of eigenstates Beugeling  _et al._
(2015b); Liu  _et al._ (2018); LeBlond  _et al._ (2019); Haque  _et al._
(2020), and in the fluctuations of off-diagonal matrix elements Khaymovich
_et al._ (2019). Since the present $k=3$ system has strong non-chaotic
features, it is not surprising that there is even stronger deviation from
Gaussian-state behavior. Perhaps unexpectedly, EEV fluctuations show a clear
power law dependence on the Hilbert space dimension, but the exponent is
markedly different from the random-matrix prediction. This demonstrates that
eigenstates of small-$k$ Bose-Hubbard chains have non-random correlations. We
examine the most straightforward ways in which physical eigenstates could be
different from ideally ergodic eigenstates, e.g., non-Gaussian coefficient
distributions and two-point correlations between the eigenstate coefficients.
We demonstrate that none of these effects explain our observed scaling. The
observed scaling exponent thus results from more subtle structures in the
eigenstate coefficients, which seems challenging to characterize or identify.

In Section II we introduce the system, notations, and the numerical procedure
for extracting EEV fluctuations $\sigma$. Section III is the heart of the
paper and reports the main results: analytically derived trace expressions and
scaling laws for $\sigma$ in the case of idealized (Gaussian) eigenstates
(III.1), and results of extensive numerics showing where these predictions
succeed and where they fail (III.2). In Section IV we examine possible reasons
for the few site systems (particularly $k=3$) deviating from the Gaussian
prediction and displaying anomalous exponents in the $D$-dependence of
$\sigma$. We conclude that the anomalous scaling is due to subtle higher-order
correlations among the eigenfunction coefficients which cannot be captured by
the 2-point correlations among eigenfunction components. Section V provides
context and concluding discussion.

The appendices present further details, data and clarification. We show the
dependence on the interaction parameter and justify the interaction value used
in the main text (Appendices A and B). We display the structure of covariance
matrices of the physical eigenstates (C). Appendices D and E contain the
derivations of the analytic results announced in the main text.

##  II Preliminaries — Systems and Observables

In this section, we introduce the system Hamiltonian (II.1), discuss how to
define the fluctuations of the eigenstate expectation values (EEVs) (II.2),
and introduce notation to be used in this paper.

###  II.1 Hamiltonians

We will investigate Bose-Hubbard systems restricted to open-boundary chains of
length $k$, with nearest neighbor hoppings and on-site interactions. The
Hamiltonian is

$H=-1/2\sum_{\langle
i,j\rangle}J_{i,j}a_{i}^{\dagger}a_{j}+\frac{U}{2}\sum_{j=1}^{k}n_{j}(n_{j}-1),$
(2)

where $\langle i,j\rangle$ denotes summation over adjacent sites ($j=i\pm 1$),
$a_{j}^{\dagger}$ and $a_{j}$ are the bosonic creation and annihilation
operators for the $j$-th site and $n_{j}=a_{j}^{\dagger}a_{j}$ is the
corresponding occupation number operator. $J_{i,j}=J_{j,i}$ is the symmetric
tunneling coefficient and $U$ is the two-particle on-site interaction
strength. We choose $J_{1,2}=1.5$ and $J_{i,j}=1$ for $i,j\geq 2$ to avoid
reflection symmetry. The particle number $N$ is conserved by the Hamiltonian.
We introduce the tuning parameter $\Lambda=UN$.

In the limits $\Lambda\to 0$ and $\Lambda\to\infty$ the system is integrable.
If $\Lambda=0$, then the free bosonic system is solvable in terms of single-
particle eigenstates, while $\Lambda\to\infty$ effectively means $J$ can be
neglected, so that Eq. (2) is diagonal already. For intermediate $\Lambda$ the
systems behave chaotically, i.e., the spectrum approximately obeys Wigner-
Dyson statistics. For the smallest chain length $k=3$, there are stronger
deviations, but nevertheless the spectral statistics near the middle of the
spectrum is close to Wigner-Dyson form (Appendix A).

The Hilbert space size is $D=\binom{N+k-1}{k-1}=\binom{N+k-1}{N}$, which
scales for fixed $k$ in the large $N$ limit as $D\sim N^{k-1}$. Therefore, for
fixed $k$ in the large $D$ limit, $N$ scales as $N\sim D^{1/(k-1)}$.

The operators we focus on are the tunnel operator $a_{2}^{\dagger}a_{1}$ from
site 1 to site 2 and the number operator $n_{1}$ on site 1. We have checked
that the overall scaling behaviors are the same for $a_{i}^{\dagger}a_{j}$
with other $i$, $j$. Note that $a_{2}^{\dagger}a_{1}$ is non-hermitian. If
hermiticity is required, we will instead investigate
$(a_{2}^{\dagger}a_{1}+a_{1}^{\dagger}a_{2})$. This change results in an
additional factor of 2 in the corresponding EEV fluctuation.

###  II.2 Defining EEV fluctuations

Figure 1:  Eigenstate expectation values $A_{\alpha\alpha}$ of the tunnel
operator $A=a_{2}^{\dagger}a_{1}$, plotted against eigenenergies, for a Bose-
Hubbard chain with a) $k=3$ sites and b) $k=10$ sites. The numbers of
particles $N$ are listed in legends. The interaction parameter is
$\Lambda\approx 2.477$ as explained in the text. Top panels show full spectra.
Bottom panels zoom into the 5th of ten equal energy intervals, as indicated by
the shorthand label “$\Delta{E}=5$”. Dotted lines are fitted linear functions.

We denote by $\left|E_{\alpha}\right\rangle$ the eigenstate of the Hamiltonian
$H$ belonging to the eigenvalue $E_{\alpha}$. If the ETH holds,
$A_{\alpha\alpha}=\left\langle
E_{\alpha}\right|A\left|E_{\alpha}\right\rangle$ is constant on small energy
windows $[E_{\alpha}-\delta,E_{\alpha}+\delta]$ around $E_{\alpha}$ for
suitable (small) $\delta$, up to finite size fluctuations. Further,
$A_{\alpha\alpha}$ coincides with the microcanonical average

$\langle A\rangle_{(E_{\alpha},\delta)}=\frac{1}{N_{E_{\alpha},\delta
E_{\alpha}}}\sum_{\beta:E_{\beta}\in[E_{\alpha}-\delta,E_{\alpha}+\delta]}A_{\beta\beta}.$
(3)

A quantitative measure of how well the EEVs coincide with the microcanonical
average is the magnitude of the finite size fluctuations, i.e., the
statistical standard deviation of $A_{\alpha\alpha}-\langle
A\rangle_{(E_{\alpha},\delta)}$,

$\left[\frac{1}{N_{\Delta E}}\sum_{E_{\alpha}\text{ in }\Delta
E}\left(A_{\alpha\alpha}-\langle
A\rangle_{(E_{\alpha},\delta)}\right)^{2}\right]^{1/2}$ (4)

over some energy window $\Delta E$, where $N_{\Delta E}$ is the number of
eigenvalues $E_{\alpha}$ in $\Delta E$. Graphically this measures the width of
the distribution of $A_{\alpha\alpha}$ values.

The above definition involves a choice of the window width $\delta$; this
choice may depend on various parameters of the physical model. To avoid larger
scale changes of $A_{\alpha\alpha}$ (e.g, as visible in Figure 1), the window
should not be too wide. The window should also not be too small, so as to
ensure that sufficiently many EEVs are within the window around $E_{\alpha}$
to get a statistically significant estimate of the fluctuations. Choosing
$\delta$ can thus be a tricky balancing act.

We can avoid these technicalities by noting that ETH implies that the quantity
$A_{\alpha\alpha}$ will vary smoothly as a function of $E_{\alpha}$, up to
finite size fluctuations. If the ETH holds in this sense, then, by definition
of smoothness, $A_{\alpha\alpha}$ should be locally linear. We divide the
energy spectrum into 10 equal-length, disjoint intervals $\Delta E$. In these
intervals the large scale change of $A_{\alpha\alpha}$ in our case is indeed
linear to an excellent approximation, as seen in the lower panels of Figure 1.
We fit linear functions $E_{\alpha}\to b+mE_{\alpha}$, on each interval to
$A_{\alpha\alpha}$. We then investigate the fluctuations around these
functions, i.e.,

$\sigma^{2}(A,\Delta E)=\frac{1}{N_{\Delta E}}\sum_{E_{\alpha}\text{ in
}\Delta E}|A_{\alpha\alpha}-b-mE_{\alpha}|^{2}.$ (5)

Since we are primarily interested in mid-spectrum eigenstates, we will show
data from the 5th, 6th and 7th energy intervals. In shorthand these will be
labeled as $\Delta{E}=5$, $\Delta{E}=6$, $\Delta{E}=7$, with $\Delta{E}$
referring to the label and not the interval width. In Figure 1, we display the
EEV’s for the full spectrum and for the $\Delta{E}=5$ interval, for two
different values of $k$. Unless indicated otherwise, we present data for an
intermediate value of the interaction parameter $\Lambda=UN$ around which the
systems are found to be significantly chaotic, namely,
$\Lambda=10^{13/33}\approx 2.477$. Appendix A provides further details on this
choice of $\Lambda$.

##  III EEV fluctuation scaling

In this section we first (III.1) consider Gaussian-random states, e.g.
eigenstates of matrices drawn from the GOE ensemble. We present expressions
for the EEV fluctuations for such idealized eigenstates. For operators of the
type $A=a_{i}^{\dagger}a_{j}$, we show that the scaling of EEV fluctuations,
$\sigma\sim{D}^{-e_{0}}$, occurs with exponent
$e_{0}=\frac{1}{2}-\frac{1}{k-1}$ in the classical limit ($N\to\infty$, fixed
$k$). For scaled operators $\bar{A}=A/N$, the exponent is $e_{0}=\frac{1}{2}$.
In comparison, in the usual thermodynamic limit ($N\to\infty$, fixed $k/N$),
the exponent is $e_{0}=\frac{1}{2}$ for local operators.

Next (III.2), we present numerically calculated scaling results for Bose-
Hubbard chains. We show that the EEV fluctuations have power-law dependence on
the Hilbert-space dimension, $\sim{D}^{-e}$. The exponent $e$ matches the
Gaussian-state prediction $e_{0}$ for larger chains, but differs substantially
for small $k$.

###  III.1 Derivation for random Gaussian states

The EEV’s calculated for random states do not have large-scale smooth
variation as a function of energy, in contrast to Figure 1. Thus, the
statistical standard deviation of EEV’s, $\sigma$, can be directly compared
with our measurement of EEV fluctuations in the physical eigenstates. The
assumption of eigenstates being effectively random has been previously used to
derive scaling properties of EEV’s in the thermodynamic limit Neuenhahn and
Marquardt (2012); Beugeling  _et al._ (2014); Khaymovich  _et al._ (2019).
Here, we provide explicit expressions for $\sigma$ in terms of trace
properties of the operator matrix, and then specialize to both thermodynamic
and classical limits. We present the important expressions here, and relegate
derivation details to Appendices D and E.

####  III.1.1 Trace expressions

Let $A$ be a $D\times D$ square matrix representing the operator of interest,
and $\left|Z\right\rangle$ be a $D$-dimensional, multivariate random state
with identically and independently distributed (i.i.d.) components $Z_{i}$,
each with mean 0. Then the statistical variance of $\langle A\rangle=\langle
Z|A|Z\rangle$ can be expressed as

$\displaystyle\sigma^{2}(\langle A\rangle)$
$\displaystyle=\left(E[Z_{1}^{4}]-3E[Z_{1}^{2}]^{2}\right)\sum_{i}A_{ii}^{2}$
$\displaystyle+E[Z_{1}^{2}]^{2}\left[\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})\right].$
(6)

Here $E[\cdot]$ represents the expectation value. The manipulations leading to
this expression can be found in Appendix D.

Eq. (III.1.1) is valid irrespective of the distribution of the i.i.d.
variables $Z_{i}$. Simplifications are achieved by specializing to the case of
GOE eigenstates, widely considered as reasonable models for the behavior of
mid-spectrum eigenstates of chaotic Hamiltonians. Eigenstates of
$D$-dimensional Gaussian matrices are uniformly distributed on the
$D$-dimensional unit sphere. In the large $D$ limit their components can be
regarded as independently normally distributed with mean 0 and variance $1/D$.
If we take $Z_{i}$ to be normally distributed with variance $1/D$, in addition
to the previous constraints, then the first term in Eq. (III.1.1) vanishes.
Thus

$\sigma^{2}(\langle
A\rangle)=\frac{1}{D^{2}}\left[\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})\right]$
(7)

for random states with Gaussian-distributed coefficients.

The second of the two traces, $\operatorname{tr}(AA^{\dagger})$, is the
squared Hilbert-Schmidt norm or the Frobenius norm of the operator. For
Hermitian $A$ the two trace terms are equal:
$\operatorname{tr}(A^{2})=\operatorname{tr}(AA^{\dagger})$. As
$\operatorname{tr}(A^{2})$ and $\operatorname{tr}(AA^{\dagger})$ are invariant
under a basis change, so is the variance Eq. (7). In contrast, Eq. (III.1.1)
is not basis-invariant, due to the first term. Expressions equivalent or
analogous to Eq. (7) have appeared previously in the literature, e.g., in
Refs. Gemmer  _et al._ (2009); Reimann (2007, 2008); Lloyd (2013).

For observables of the form $A=a_{j}^{\dagger}a_{i}$, the trace expression is
shown in Appendix E to be given by

$\displaystyle\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})$
$\displaystyle=\frac{N(N+k)}{k(k+1)}D$ $\displaystyle i\neq j,$ (8a)
$\displaystyle\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})$
$\displaystyle=\frac{2N(2N+k-1)}{k(k+1)}D$ $\displaystyle i=j.$ (8b)

In the usual thermodynamic limit, the fraction multiplying $D$ has $O(D^{0})$
scaling in either case. Thus one obtains $\sigma^{2}\sim\frac{D}{D^{2}}$,
i.e., the usual $\sigma{\sim}D^{-1/2}$ scaling for EEV fluctuations in the
thermodynamic limit.

####  III.1.2 Scaling in the classical limit

We now consider the classical limit, $k\ll N$. If $A$ is a linear combination
of terms like $a_{j}^{\dagger}a_{i}$, the trace expression scales as

$\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})\sim N^{2}D.$ (9)

Since $N\sim D^{1/(k-1)}$ in the classical limit, the variance scales as

$\sigma^{2}(\langle A\rangle)\sim\frac{D\cdot
D^{2/(k-1)}}{D^{2}}=D^{-2e_{0}},$ (10)

where

$e_{0}=e_{0}(k)=\frac{1}{2}-\frac{1}{(k-1)}$ (11)

is the scaling exponent announced previously. For $k\gg 1$, but still
$k\ll{N}$, the second term becomes negligible and we obtain
$\sigma\sim{D}^{-1/2}$ scaling similar to the thermodynamic limit.

For numbers of sites that are not too large, the EEV scaling of two-point
operators (of the type $a_{j}^{\dagger}a_{i}$ or their linear combinations) is
different for the classical limit compared to the thermodynamic limit.
Mathematically, this difference can be attributed to the operators $A$ scaling
with $N$. If we rescale $A$ to $\bar{A}=A/N$, then the traces (9) scale as $D$
rather than as $N^{2}D$, so that the variance scales for all $k$ as

$\sigma^{2}(\langle\bar{A}\rangle)\sim\frac{D}{D^{2}}=D^{-1}.$ (12)

Summarizing: in the classical limit, the EEV fluctuation scaling is
${\sim}D^{-1/2}$ for scaled operators $\bar{A}=A/N$ for all $k$, and also for
un-scaled operators $A$ in the $k\gg 1$ limit. This is the same exponent
$e_{0}=\frac{1}{2}$ familiar from the thermodynamic limit Beugeling  _et al._
(2014); Khaymovich  _et al._ (2019). However, for moderate $k$ and for the
operator $A$, the scaling is according to the exponent
$e_{0}=\frac{1}{2}-\frac{1}{(k-1)}$ of Eq. (11)

###  III.2 Numerical results: Bose-Hubbard eigenstates

Figure 2:  EEV fluctuations $\sigma$ of the operator $A=a_{2}^{\dagger}a_{1}$
against Hilbert space size $D$ for various chain lengths $k$. Data shown for
eigenstates in energy windows $\Delta E=$ 5 (blue circles), 6 (red squares), 7
(green inverted triangles). The $\sigma$ vs $D$ data sequences are arranged
reasonably linearly in all cases in the log-log plots, suggesting
$\sigma\sim{D}^{-e}$ behavior. The slopes of fitted lines (i.e., numerical
estimates of $-e$) are given in the legends. The Gaussian predictions $e_{0}$
for the exponents are $0$, $\frac{1}{6}\approx 0.1667$, $0.25$, and
$\frac{7}{18}\approx 0.3889$ respectively for $k=$ 3, 4, 5, and 10.

In Figure 2 we show the fluctuations $\sigma$ of the EEVs for different energy
windows near the middle of the spectrum, plotted against $D$. Each panel shows
a different (fixed) number of sites $k$; in each case the classical limit is
approached by increasing $N$. Generally, the sequences follow clear power-law
dependencies, $\sigma\sim{}D^{-e}$. The power-law behavior sets in at
relatively small values of $D$ already.

It is clear from the $k=3$ data, panel (a), that the exponent $e$ does not
match the value predicted for Gaussian-random states, Eq. (11), which is
$e_{0}=0$ for $k=3$. The EEV fluctuation for the system eigenstates
_increases_ with a positive exponent ($e<0$) instead of being flat as a
function of $D$. Similarly, for the $4$-site chain the exponent $e$ is seen to
be slightly negative — $\sigma$ increases slowly with system size — whereas
the predicted value is $e_{0}=+1/6$.

Figure 3:  The exponent $e$ versus the chain length $k$ for Bose-Hubbard
eigenstates in different energy windows, and for Gaussian-random states. The
pink dashed curve is the predicted formula for exponents, $e_{0}$, which tends
to $1/2$ (solid horizontal line) for large $k$. Inset to left panel shows the
error bars for the estimation of $e$ from $\Delta{E}=5$ data. The error bars
are omitted elsewhere and will be omitted in later figures.

The calculations rely on full numerical diagonalization, and hence are limited
by the Hilbert space size $D$. Our limit was $D\lesssim 10^{5}$. For each $k$,
we increased the particle number $N$ as far as possible such that $D$ did not
exceed $100,000$. For small $k$, this provides a satisfactory number of
available $N$ values, and extracting the exponent $e$ from a fit to
$\sigma\sim{}D^{-e}$ is quite reliable. For large $k$, only a few $N$ values
are available. For the largest lattice ($k=15$), only three data points
($N=4$, $N=5$, and $N=6$) were used. This means a large uncertainty in the
estimation of $e$ (Figure 3 inset). It also means that the regime $N\gg k$ is
not reached.

In Figure 3 we present the exponents $e$ extracted from the numerical data. In
addition to the exponents for the operator $a_{2}^{\dagger}a_{1}$
(corresponding to Figure 2), we also show the exponents for the operator
$n_{1}$. For small $k$, the numerically observed exponents $e$ fall
significantly below the Gaussian-random case, for both operators. For larger
$k$ values, the Bose-Hubbard systems show EEV fluctuations closer to the
Gaussian-random case, at least for $\Delta{E}=5,6$. (The $\Delta{E}=7$ window
shows larger deviation, presumably becuase it is closer to the edges of the
spectrum.) We interpret this as a signature of the large-$k$ Bose-Hubbard
systems being more chaotic, so that mid-spectrum eigenstates are better
approximated by Gaussian-random states. The deviation for small $k$ represents
the ‘mixed’ (chaotic+regular) nature of the few-site Bose-Hubbard
Hamiltonians.

Figure 3 also shows numerically calculated exponents for EEV fluctuations in
Gaussian-random states (pink triangles), and compares with the $k\ll{N}$
prediction, Eq. (11) (pink dashed curve). The agreement is good for all $k$
and excellent for small $k$. At larger $k$, computer memory limitations
prevent our calculations from reaching particle numbers $N\gg{k}$. This
explains the (minor) deviation of the numerical exponents from the $N\gg{k}$
prediction.

Figure 4:  Similar data as in Figures 2 and 3 but with scaled operators
$\bar{A}$. The operators are scaled by factors $\sim{N}$; the precise factors
are explained in the text.

One can view the same effects through the fluctuations of the scaled
($\bar{A}$) operators $a_{2}^{\dagger}a_{1}/N$ and $n_{1}/N$. For these
operators, the predicted exponent is $1/2$ for all $k$. We present the
numerical exponents for such operators in Figure 4; however we scale by
factors slightly different from $N$. The prediction $e_{0}=1/2$ was obtained
in the previous subsection by assuming $N\gg k$. In particular in the trace
expressions of Eq. (8), this led to $N(N+k)\approx{N^{2}}$ and
$2N(2N+k-1)\approx{4N^{2}}$. In our numerical calculations, for larger $k$ we
do not have access to $N$ values in this regime. Therefore we scale the
operators by the factor $\sqrt{N(N+k)}$ for $A=a_{i}^{\dagger}a_{j}$ with
$i\neq{j}$ and by the factor $\sqrt{N(2N+k-1)}$ for $i=j$. With this
modification, the numerically calculated exponents using Gaussian-random
states (pink triangles) do not deviate systematically from $1/2$ at large $k$,
even though the $N\gg{k}$ regime is not reached. The observed physical
exponents (for Bose-Hubbard eigenstates) are significantly different from the
predicted $e_{0}=1/2$ for small $k$, but approach this value as $k$ is
increased.

Summarizing our numerical findings: power-law dependence of EEV fluctuations
on $D$ is seen for all $k$. The exponent at larger $k$ (more fully chaotic
systems) matches well the Gaussian prediction. For smaller $k$ (mixed
systems), there are significant deviations from the Gaussian prediction.
Perhaps the most remarkable aspect of these numeric results is that, for small
$k$, the departure of the eigenstates from Gaussian-random or ergodic behavior
does not destroy the power-law dependence of EEV fluctuations with $D$ but
changes the exponent substantially. In the following section we will examine
possible explanations for this phenomenon.

##  IV Non-reasons for anomalous scaling

The mismatch with predicted scaling for small $k$ must be due to the
eigenstates of few-site Bose-Hubbard systems deviating from idealized
Gaussian-random states, which we refer to as ‘ergodic’ states. In this
section, we examine various types of deviation from the idealized case, and
rule out several plausible explanations for the anomalous scaling.

One could imagine that the eigenstates effectively occupy a smaller fraction
of the Hilbert space than a random state, and that this fraction has sublinear
scaling with $D$. This can be quantified through analysis of the participation
ratio or multifractal dimensions of the eigenstates. Therefore, in Subsection
IV.1, we discuss the participation ratio $P_{\alpha}$. We show how
$P_{\alpha}\sim D^{1}$ scaling implies the expected EEV fluctuation scalings
that we have derived previously in III.1. We also show that Hilbert space
occupancy is not the explanation for our anomalous EEV scaling exponents,
because eigenstates of the $k=3$ system have $P_{\alpha}\sim D^{1}$ scaling.

One could also imagine that the anomaly of scaling exponents is due to the
eigenstate components not being identically distributed. In IV.2 we present
data showing that this is not the reason for the anomalous scaling of EEV
fluctuations.

These results show that the non-ergodic scaling must be due to correlations
present in the eigenstates. In IV.3 we examine the simplest (and most
prominent) type of correlations between eigenstate coefficients, namely, those
captured in the covariance matrix. Perhaps surprisingly, we find that these
correlations do not explain the anomalous scaling either. The anomalous
scaling exponent for $k=3$ EEV fluctuations is thus caused by more subtle
(higher-order) correlations.

###  IV.1 Participation Ratios

For this Subsection it will be convenient to consider our operator $A$ to be
Hermitian. We expand the eigenstates $\left|E_{\alpha}\right\rangle$ of the
Hamiltonian in the eigenstate basis $\left|\phi_{\gamma}\right\rangle$ of the
operator $A$

$\left|E_{\alpha}\right\rangle=\sum_{\gamma}c_{\gamma}^{(\alpha)}\left|\phi_{\gamma}\right\rangle,$
(13)

where
$c_{\gamma}^{(\alpha)}=\left\langle\phi_{\gamma}|E_{\alpha}\right\rangle$. If
we denote the eigenvalues of $A$ as $a_{\gamma}$, then the EEVs can be written
as

$A_{\alpha\alpha}=\sum_{\gamma=1}^{D}|c_{\gamma}^{(\alpha)}|^{2}a_{\gamma}.$
(14)

We regard the coefficients $c_{\gamma}^{(\alpha)}$ to be random variables,
with each eigenstate index $(\alpha)$ denoting a different sample from the
same underlying distribution. The distribution is assumed to be the same for
every $\gamma$. As usual, this framework will not capture large-scale
dependences of the EEVs $A_{\alpha\alpha}$ on the energies $E_{\alpha}$. This
is acceptable because we are interested in the fluctuations only.

The variance of the EEVs is

$\operatorname{var}\left(\sum_{\gamma=1}^{D}|c_{\gamma}|^{2}a_{\gamma}\right)=\sum_{\gamma=1}^{D}\operatorname{var}\left(|c_{\gamma}|^{2}\right)a_{\gamma}^{2}\\\
=\operatorname{var}\left(|c_{\gamma}|^{2}\right)\operatorname{tr}(A^{2}).$
(15)

The variance of $|c_{\gamma}|^{2}$ can be written as

$\operatorname{var}\left(|c_{\gamma}|^{2}\right)=E\left[|c_{\gamma}|^{4}\right]-\left(E\left[|c_{\gamma}|^{2}\right]\right)^{2}=\frac{1}{DP}-\frac{1}{D^{2}}$
(16)

where we have used the definition of the participation ratio to be

$P_{\alpha}=\left(\sum_{\gamma=1}^{D}|c_{\gamma}^{(\alpha)}|^{4}\right)^{-1}=\left(\sum_{\gamma=1}^{D}D\times
E\left[|c_{\gamma}^{(\alpha)}|^{4}\right]\right)^{-1}\ .$ (17)

We thus have the prediction for the EEV variance

$\sigma^{2}=\left(\frac{1}{DP}-\frac{1}{D^{2}}\right)\operatorname{tr}(A^{2})=\left(\frac{1}{DP}-\frac{1}{D^{2}}\right)\operatorname{tr}(AA^{\dagger})$
(18)

for Hermitian operators. For _Gaussian_ states, $P=D/3$, so that this
expression for $\sigma^{2}$ reduces to Eq. (7). More generally, as long as the
participation ratio scales linearly with $D$, the factor in brackets $\sim
1/D^{2}$, so that we obtain the same scaling as for Gaussian states. For non-
ergodic states $P\sim D^{K}$ with $K<1$, the first term in brackets would
dominate and one would obtain different scaling.

Figure 5:  Participation ratio in the basis of eigenstates of the operator
$a_{2}^{\dagger}a_{1}+a_{1}^{\dagger}a_{2}$. Horizontal lines indicate the
Gaussian expectation ($=\frac{1}{3}$). Top panels: Scaled participation ratio
$P_{\alpha}/D$ of energy eigenstates, versus corresponding energy eigenvalues.
Lower panels: average scaled participation ratio for different energy
intervals plotted against Hilbert space dimension $D$, for fixed chain length
$k$ and increasing particle number $N$.

In Figure 5 we show the behavior of the participation ratio in the basis of
eigenstates of $(a_{2}^{\dagger}a_{1}+a_{1}^{\dagger}a_{2})$. The mid-spectrum
$P_{\alpha}$ is close to the Gaussian expectation for highly chaotic (larger
$k$) systems. For $k=3$ the deviation from Gaussianity ($P=D/3$) is strong.
However, in both cases the scaling of $P_{\alpha}$ with Hilbert space
dimension $D$ is very linear (lower panels).

Thus, Hilbert space occupancy does not explain the observed anomalous exponent
of EEV fluctuations.

###  IV.2 Non-identical distribution of eigenstate coefficients

The analysis in Subsections IV.1 and III.1 are based on eigenstate
coefficients being _identically distributed_ and _independent_. Could it be
that the failure to capture the EEV fluctuation scaling at small $k$ results
from the eigenstate components having non-identical distributions?

To check, we drop the condition that the eigenstates have identically
distributed components and only assume that the components are independent. In
Figure 6 we show the EEV fluctuations obtained from an estimation of the
underlying distributions of the eigenstate components of the Bose-Hubbard
systems. We sampled uniformly from components of all eigenstates of the Bose-
Hubbard model in a specific energy interval. Effectively, for each component
of the state, we randomly picked an eigenstate in the energy interval and used
the entry of the corresponding component. The EEV fluctuations are then
calculated from such sampled states. The resulting data are marked
‘independent’ in Figure 6. The results match well with the EEV fluctuations
obtained from Gaussian states and do not match the EEV fluctuations of the
actual physical systems for small $k$.

We conclude that that any model for the wavefunctions with independent
components does not explain the observed anomalous scaling at small $k$. In
other words, the root cause of the phenomenon is not the eigenstate
coefficients having non-identical or non-Gaussian distributions, or
insufficient occupancy of the Hilbert space. Rather, we have traced the cause
to the fact that the eigenstate coefficients are not independent.

Figure 6:  The operator in all plots is $a_{2}^{\dagger}a_{1}$. Blue dots:
eigenstates of physical Bose-Hubbard system. Red squares: Gaussian states with
i.i.d. components. Green down triangles: vectors with independent but non-
identically distributed components, each component sampled from system
eigenstates. Purple up triangles: multivariate Gaussian states with covariance
matrix estimated from system eigenstates. a)-b) EEV fluctuations. The dashed
purple line is the prediction by Eq. (20). c) The exponent $e$, such that
$\sigma\sim{D}^{-e}$, versus the number of sites $k$ for the same
distributions as in a)-b).

###  IV.3 Eigenstate Correlations

Continuing our effort to identify what feature of small-$k$ eigenstates is
responsible for the anomalous ETH scaling, we relax the constraints on the
model states even further. Now, we assume that eigenstates are drawn
independently from a multivariate distribution $Z$, but the $D$ components of
$Z$ are allowed to be correlated. Given a multivariate distribution $Z$, the
statistical correlations between components are quantified by the covariance
matrix $\Sigma$. The covariance matrix in the case of eigenstates can be
estimated by regarding the eigenstates within a specific energy interval as
different samples of $Z$.

It is reasonable to assume that the mean of all components of $Z$ is zero, as
components of mid-spectrum eigenstates generally have zero mean. Then an
estimate of $\Sigma$ is given by

$\Sigma=N_{\Delta E}^{-1}\sum_{E_{\alpha}\text{ in }\Delta
E}\left|E_{\alpha}\right\rangle\left\langle E_{\alpha}\right|,$ (19)

where $E_{\alpha}$ and $\left|E_{\alpha}\right\rangle$ denote eigenvalues and
eigenstates, respectively, and $N_{\Delta{E}}$ is the number of eigenstates in
the energy window $\Delta{E}$. Eq. (19) follows directly from the definition
of the sample covariance matrix of $Z$, where the samples are the eigenstates
$\left|E_{\alpha}\right\rangle$.

In Appendix C, Figure 9, we have provided visualizations of estimated
covariance matrices for $k=3$ and $k=6$. The $k=3$ case is seen to have
significant off-diagonal elements (arranged in intriguing patterns),
indicating non-negligible correlations between the eigenstate components.

Once we have constructed the correlation matrix $\Sigma$ fom the actual
eigenstates, we can sample $D$-component vectors whose components are
distributed _and_ correlated according to $\Sigma$. In Figure 6 the EEV
fluctuations obtained from such sampled states are marked as ‘correlated’. The
values of the fluctuations thus obtained are larger than those obtained from
the independent-component random states, and more comparable to the
fluctuations obtained from the physical eigenstates. In the chaotic cases
(large $k$), all of these cases have the same scaling. However, in the $k=3$
case, the scaling exponent is close to the ergodic (Gaussian) case and does
not reproduce the anomalous scaling at all. This is seen in panel a) of Figure
6, and also in panel c) where the fitted exponents are plotted. The fitted
exponent is slightly off the Gaussian value for small $k$, but far from the
anomalous values of the physical system.

These results show that the deviation of the physical system from expected
‘ergodic’ behavior is only partially captured by the two-point correlations
between the eigenstate components. This suggests that the small-$k$
eigenstates deviate from randomness in some more drastic manner, which does
not seem easy to quantify.

In addition to the direct numerical verification discussed above, one can use
analytic results to argue (non-rigorously) that inclusion of reasonable two-
point correlations in the model of random states should not change the EEV
fluctuation scaling exponent. Given a multivariate Gaussian distribution with
mean zero and covariance matrix $\Sigma$, we can show (details in Appendix D)
that the variance of EEV fluctuations is given by

$\operatorname{var}(\langle
A\rangle)=\frac{1}{D^{2}}\left[\operatorname{tr}((D\Sigma
A)^{2})+\operatorname{tr}(D\Sigma AD\Sigma A^{t})\right].$ (20)

This is Eq. (7) with the change $A\to D\Sigma A$. The variance of components
of the wavefunction is not fixed by normalization of the eigenstates any more.
But it is reasonable to assume that it still scales as $\sim 1/D$, and we have
checked numerically that this scaling holds for the mid-spectrum eigenstates
of all our physical systems, including $k=3$. Since the variances of the
wavefunction components are the diagonal entries of $\Sigma$, the diagonals of
$D\Sigma$ scale (at most) as constant in $D$. By the Cauchy-Schwartz theorem
the off-diagonal terms of $\Sigma$ are bounded by the diagonal, so every
component of $D\Sigma$ is (at most) constant in $D$.

Without making assumptions about the detailed structure of $\Sigma$, we cannot
derive rigorously the scaling of the traces in Eq. (20), which was possible
for Eq. (8) or (9). However, since $D\Sigma$ is elementwise at most
$\sim{D^{0}}$, and assuming $\Sigma$ is not too exotic, one can argue that the
derivation of Eq. (E) should hold for this case as well. In other words, for
‘reasonable’ $\Sigma$, one expects the same scaling as in the case of
independent Gaussian eigenstates. This is consistent with Figure 6, where the
$\Sigma$ is estimated numerically from the physical eigenstates.

##  V Discussion and Context

Motivated by questions relating quantum dynamics to statistical mechanics, we
have undertaken a study of the eigenstate thermalization hypothesis in the
scaling sense, but considering increasing Hilbert space dimensions along the
(semi)classical limit rather than the usual thermodynamic limit. This has led
to a characterization of the distinctive properties of few-site Bose-Hubbard
systems in terms of anomalous scaling exponents.

Summary of analytic results. For GOE eigenstates, i.e., states whose
components can be approximated by Gaussian-distributed random variables, we
have used trace expressions for the EEV fluctuation $\sigma$, Eq. (7). For
operators of the type $A=a_{j}^{\dagger}a_{i}$, the trace operators can be
expressed as Eq. (8). Based on these main expressions, we are able to predict
ideal scaling behaviors of EEV fluctuations in the classical limit, for both
un-scaled operators of the type $A$ and scaled operators $\bar{A}=A/N$. Of
course, the usual ETH scaling of the thermodynamic limit also follows from
these expressions.

In the classical limit $N\gg{k}$, the EEV fluctuations are found for such
idealized eigenstates to behave as $\sigma\sim{D}^{-e_{0}}$, with
$e_{0}=\frac{1}{2}-\frac{1}{(k-1)}$ for un-scaled $A$ operators and
$e_{0}=1/2$ for scaled $\bar{A}$ operators.

In addition, we have presented expressions for $\sigma$ for a number of
related cases, e.g., for i.i.d. distributed eigenstate components with the
distribution not assumed to be Gaussian, Eq. (III.1.1), in terms of the
participation ratio, Eq. (18), and for the more general case where the
eigenstate components are allowed to be correlated according to a covariance
matrix, Eq. (20).

Summary of numerical results. We have explored the scaling exponent for
various lattice lengths $k$, increasing the boson number $N$ with fixed $k$ to
approach the classical limit. At larger $k$, the exponent matches the random-
eigenstate prediction. At small $k$, the fluctuation appears to have power-law
dependence $\sigma\sim{D}^{-e}$ on the Hilbert space dimension, i.e., $e$ is
well-defined, but the value of $e$ differs markedly from the ergodic
prediction. Through a series of additional numerical tests, we have shown that
this anomalous scaling is not explained by 2-point correlations between
eigenstate components.

The small-size Bose-Hubbard systems thus have mid-spectrum eigenstates which
violate the usual randomness approximation in some subtle higher-order manner.

Deviation from Ergodicity. The deviation of quantum many-body systems from
ergodicity has been the subject of interest from multiple viewpoints in recent
years. Other than the strong violations of ETH/ergodicity due to integrability
or many-body localization, more subtle departures have also been addressed or
observed. In many-body systems that are nominally chaotic, mid-spectrum states
are largely well-modeled by random states, but small or subleading deviations
have been observed in various properties Beugeling  _et al._ (2018); Luitz and
Bar Lev (2016); Bäcker  _et al._ (2019); Luitz  _et al._ (2020); Beugeling
_et al._ (2015b); Liu  _et al._ (2018); LeBlond  _et al._ (2019); Haque  _et
al._ (2020); Khaymovich  _et al._ (2019). However, _scaling_ properties in
these systems generally follow random-state predictions. In the small-$k$
Bose-Hubbard systems, we have shown a striking exception: a system which is
not integrable or many-body localized, but nevertheless violates the usual
scaling behavior expected in chaotic systems.

The three-site and four-site Bose-Hubbard systems are widely known to be
imperfectly chaotic, in particular because their classical phase space has
been explored and is known to have both chaotic and regular regions Mossmann
and Jung (2006); Hiller  _et al._ (2009); Viscondi and Furuya (2011); Arwas
_et al._ (2014); Han and Wu (2016); Bürkle and Anglin (2019); Bychek  _et al._
(2020); Rautenberg and Gärttner (2020). With the current interest in chaos in
many-body eigenstates, it is important to characterize such deviations from
chaoticity. In this work, ETH scaling (i.e., the scaling of EEV fluctuations)
has turned out to be fruitful approach to characterize these special systems.

Open questions and ideas arising from this work. The present work opens up a
number of new questions deserving investigation:

(1) We have found that the small-$k$ Bose-Hubbard systems display EEV
fluctuations scaling with exponents that appear numerically well-defined but
very clearly different from the random-state prediction. An analytic
explanation for these observed new exponents is currently not available, and
remains an open question. The anomalous scaling is related to the insufficient
chaoticity of few-site systems, which is closely connected to the mixed phase
space of the corresponding classical system. Hence a tempting conjecture is
that some property measuring the degree of chaos in the classical limit might
explain the exponents.

(2) Our analytic results have focused on essentially infinite-temperature
states. It would be interesting to develop trace expressions for finite
temperatures. This is likely not possible to do in complete generality without
making assumptions on the system Hamiltonian, but perhaps some results can be
derived with minimal assumptions, such as locality of the Hamiltonian.

(3) For the thermodynamic limit, the EEV fluctuations of _integrable_ systems
generally show power-law scaling in $N$ or $k$ Ziraldo and Santoro (2013);
Ikeda  _et al._ (2013); Beugeling  _et al._ (2014); Alba (2015); Nandy  _et
al._ (2016); Magán (2016); Haque and McClarty (2019); Mierzejewski and Vidmar
(2020), which translates into $\sim\ln{D}$ scaling. For the classical limit,
however, power-law scaling in $N$ would mean power-law scaling in $D$ as well.
Our numerics (Appendix B) shows that fitting $\sigma\sim{D}^{-e}$ in the
integrable regime yields $e=-1/(k-1)$ for $A$-type operators and $e=0$ for
$\bar{A}$-type operators. A detailed understanding and explanation of these
exponents remains a task for future research.

(4) Bose-Hubbard systems are, of course, not the only quantum systems with a
classical limit. It remains to be discovered how generic our findings are.
Comparing chaos-related properties between quantum systems and the
corresponding classical systems has also been of interest in few-spin systems
Feingold and Peres (1983); Feingold  _et al._ (1984); Robb and Reichl (1998);
Emerson and Ballentine (2001); Ballentine (2004); Fan  _et al._ (2017);
Pappalardi  _et al._ (2018); Piga  _et al._ (2019); Lerose and Pappalardi
(2020); Xu  _et al._ (2020); Pappalardi  _et al._ (2020); Mondal  _et al._
(2020) and in spin-boson systems Bastarrachea-Magnani  _et al._ (2014, 2015);
Ray  _et al._ (2016b); Bastarrachea-Magnani  _et al._ (2016); Lewis-Swan  _et
al._ (2019); Chávez-Carlos  _et al._ (2019); Lerose and Pappalardi (2020); Xu
_et al._ (2020); Wang and Robnik (2020); Villaseñor  _et al._ (2020);
Pilatowsky-Cameo  _et al._ (2020); Sinha and Sinha (2020). (The classical
system is obtained in the limit of large spin quantum number.) Studying the
behavior of EEV fluctuations in such systems when approaching the classical
limit would provide interesting characterizations of ergodicity, e.g., of how
well randomness approximations work.

(5) As part of our effort to address the anomalous scaling at small $k$, we
have briefly examined the covariance matrix of eigenstates, treating each
eigenstate as a sample drawn from the distribution of eigenstates, according
to Eq. (19). Studying the thus-defined covariance matrix might be fruitful for
various quantum systems, as the departure of this matrix from the identity
matrix tells us how different the eigenstates are from infinite-temperature
states. A further significance of this covariance matrix is that the same
object is the microcanonical density matrix, and hence its structure should
provide insights into the connection betwen quantum eigenstates and
thermodynamics.

###### Acknowledgements.

MH thanks P. Ribeiro and R. Mondaini for enlightening discussions. We thank J.
D. Urbina and K. Richter for useful correspondence. Some computations were
performed on the Irish Center for High-End Computing. GN is supported by the
Irish Research Council Government of Ireland Postgraduate Scholarship Scheme
(GOIPG/2019/58).

##  Appendix A Level spacing statistics

Figure 7:  Top: Level ratios $\tilde{r}$ averaged over the whole spectrum, as
a function of interaction strength $\Lambda=UN$, shown in the range
$\Lambda\in(0.1,100)$. The horizontal lines are $\langle\tilde{r}\rangle_{\rm
integr}$ and $\langle\tilde{r}\rangle_{\rm GOE}$. The vertical line indicates
$\Lambda=10^{13/33}\approx 2.477$. Bottom: Mean level ratio against ten evenly
spaced energy windows $\Delta E$, labeled 1 through 10 and plotted against
these labels, for $\Lambda\approx 2.477$, various system sizes. The match to
the GOE value is better away from spectral edges.

The distribution of energy eigenvalue spacings,
$s_{\alpha}=E_{\alpha+1}-E_{\alpha}$ for ordered $E_{\alpha}<E_{\alpha+1}$, is
often used as an indicator for quantum chaos. Measuring this distribution
involves unfolding the spectrum. This is bypassed by investigating the
distribution of spacing ratios Oganesyan and Huse (2007); Atas  _et al._
(2013)

$r_{\alpha}=\frac{s_{\alpha+1}}{s_{\alpha}}\quad\text{and}\quad\tilde{r}_{\alpha}=\min\left(r_{\alpha},\frac{1}{r_{\alpha}}\right).$

The latter has the advantage that it is bounded. Quantum Hamiltonians are
considered chaotic when the distributions follow that of the relevant Wigner-
Dyson ensemble, which in our case is the GOE ensemble. The mean $\tilde{r}$
for GOE is $\langle\tilde{r}\rangle_{\rm GOE}\approx 4-2\sqrt{3}\approx
0.53590$. For an integrable Hamiltonian, the level spacing distribution is
Poissonian (exponential); this leads to $\langle\tilde{r}\rangle_{\rm
integr}=2\log 2-1\approx 0.38629$.

Figure 7 (top) shows that, with increasing interaction strength $\Lambda=UN$,
the Bose-Hubbard systems turn from integrable to chaotic and back to
integrable. Some anomalies are visible for short chains (small $k$). First,
the distributions of the level ratios for very small $\Lambda$ or very large
$\Lambda$ match neither GOE nor Poisson statistics. Also, for small $k$, in
the most chaotic regime, the peak of $\langle\tilde{r}\rangle$ is still
noticeably lower than the GOE value $\approx 0.534$. As $k$ increases, the
peak gets closer and closer to the GOE value, consistent with the intuition
that larger-$k$ chains are more chaotic.

Our grid of $\Lambda$ values contains 100 values spaced logarithmically from
0.1 to 100. Among these values, $\Lambda=10^{13/33}\approx 2.477$ is found to
be close to the location of the chaotic peak, for all $k$. Therefore, in this
work we have shown data for this value of the interaction parameter.

Of course, we do not expect that the full spectrum obeys GOE statistics — the
spectral edges are expected to deviate. This can be observed in Eq. (7)
(bottom panel), where we plot $\langle\tilde{r}\rangle$ against different
energy windows.

##  Appendix B Dependence of EEV Scaling on interaction strength $\Lambda$

Figure 8:  Exponent $e$ obtained from $\sigma\sim{D}^{-e}$ fits, as a function
of the interaction parameter $\Lambda=UN$. Results are shown for both un-
scaled operator $A=a_{2}^{\dagger}a_{1}$ and its scaled version $\bar{A}$,
with $\bar{A}=A/N$ for the $k=3$ system (top) and $\bar{A}=A/\sqrt{N(N+k)}$
for the $k=6$ system (bottom).

In the main text we focused on the value of the interaction parameter
$\Lambda=UN$ at which our systems are most chaotic. Here, we show how changing
$\Lambda$ affects the scaling of EEV fluctuations. Fitting the fluctuation
data to $\sigma\sim D^{-e}$, we extract exponents $e$ which we show in Figure
8 as a function of $\Lambda$. We compare the smallest Bose-Hubbard chain with
$k=3$ sites with a larger, thus more chaotic chain with $k=6$.

We show data both for the operator $A=a_{2}^{\dagger}a_{1}$ and the scaled
operator $\bar{A}\approx a_{2}^{\dagger}a_{1}/N$. For $k=3$, we show results
for $\bar{A}=a_{2}^{\dagger}a_{1}/N$. (As the data for small $k$ extends into
the $N\gg{k}$ regime, it is reasonable to approximate $\sqrt{N(N+k)}\approx
N$.) For larger $k$ values, the regime $N\gg{k}$ is difficult to reach
numerically. Therefore it is appropriate to scale with the factor
$\sqrt{N(N+k)}$ instead of the factor $N$, as the approximation
$\sqrt{N(N+k)}\approx N$ made for Eq. (8) or Eq. (43) may not be accurate.

Since $N\sim D^{1/(k-1)}$, the exponent $e$ for $\sigma(\bar{A})$ is shifted
upwards compared to that for $\sigma(A)$ by $1/(k-1)$. This is seen in Figure
8 for both $k=3$ (shift of $\frac{1}{3-1}=0.5$) and for $k=6$ (shift of
$\frac{1}{6-1}=0.2$).

At intermediate values of $\Lambda$, when the system is (partially) chaotic,
the value of $e$ is near that discussed in the main text. If we had fully
chaotic (Gaussian) eigenstates, we would expect $e=0.5$ for the scaled
operators $\bar{A}$ and $e=0.5-0.5=0$ ($k=3$) or $e=0.5-0.2=0.3$ ($k=6$) for
the un-scaled operators $A$. Compared to the Gaussian prediction, the observed
exponents are significantly smaller for the $k=3$ system, and somewhat smaller
for the $k=6$ system, as discussed in detail in the main text.

For small or large $\Lambda$, the system is (near-)integrable. In these
limits, the fluctuations for scaled operators $\bar{A}$ do not have power-law
dependence on $D$, and the fit to $D^{-e}$ results in $e=0$. Thus eigenstate
fluctuation scaling of _scaled_ operators $\bar{A}=A/N$ in the classical limit
is analogous to that of _local_ operators in the thermodynamic limit: $e\to 0$
in the integrable limits and $e\approx 0.5$ in the chaotic regime Beugeling
_et al._ (2014). For local operators in the thermodynamic limit, the
fluctuations have power-law dependence on system size Ziraldo and Santoro
(2013); Ikeda  _et al._ (2013); Beugeling  _et al._ (2014); Alba (2015); Nandy
_et al._ (2016); Magán (2016); Haque and McClarty (2019); Mierzejewski and
Vidmar (2020), i.e., logarithmic dependence on $D$. A detailed study of the
(near-)integrable models in the classical limit remains an interesting task
for future studies.

Because $e$ for $\bar{A}$ settles to zero for $\Lambda\to 0$ and
$\Lambda\to\infty$, the exponent for the unscaled operator
$A=a_{2}^{\dagger}a_{1}$ settles to $-1/(k-1)$ in these limits, i.e., to
$-0.5$ for the trimer and to $-0.2$ for $k=6$. This is clearly seen in Figure
8.

##  Appendix C Estimated covariance matrices from physical eigenstates

Figure 9:  Scaled estimated covariance matrix $D\Sigma$ of mid-spectrum
eigenstates (energy window $\Delta E=5$) for Bose-Hubbard chains with $k$
sites and $N$ particles. The estimate $\Sigma$ is defined in Eq. (19), The
absolute values of matrix entries, $|D\cdot\Sigma_{ij}|$, corresponding to
basis $\mathcal{B}$ described in text, are shown. The arrows indicate the
ordering of the basis, e.g. in the top left is the entry
$|D\cdot\Sigma_{11}|$. $N$ had to be chosen significantly smaller than in the
rest of the paper to visualize the patterns without zooming. The patterns are
stable for increasing $N$.

In Subection IV.3 we discussed results incorporating the 2-point correlations
between eigenstate components, as embodied in the covariance matrix. The
covariance matrix for physical mid-spectrum eigenstates is calculated
according to Eq. (19), i.e., as an equal-weight mixture of the density
matrices $\left|E_{\alpha}\right\rangle\left\langle E_{\alpha}\right|$
corresponding to each eigenstate in the energy window. Here we show the
structure of the covariance matrices obtained in this way, i.e., we visualize
the 2-point correlations present in the physical eigenstates.

To visualize the structure of covariance matrices one has to fix a basis of
the Hilbert space. In Eq. (9) we show covariance matrices with respect to the
defining basis $\mathcal{B}$. The elements of $\mathcal{B}$ are the mutual
eigenstates of all number operators $n_{j}=a_{j}^{\dagger}a_{j}$, denoted by
$\left|n_{1},\dots,n_{k}\right\rangle$ where $n_{1}+\dots+n_{k}=N$. For two
states $\left|\psi\right\rangle=\left|n_{1},\dots,n_{k}\right\rangle$ and
$\left|\phi\right\rangle=\left|m_{1},\dots,m_{k}\right\rangle$ we say that
$\left|\psi\right\rangle<\left|\phi\right\rangle$ if and only if the states
interpreted as $(N+1)$-adic numbers satisfy $(\psi)_{N+1}<(\phi)_{N+1}$, i.e.

$\left|\psi\right\rangle<\left|\phi\right\rangle\iff\sum_{j=1}^{k}n_{j}(N+1)^{j}<\sum_{j=1}^{k}m_{j}(N+1)^{j}.$
(21)

The basis states in $\mathcal{B}$ are in descendent order with respect to the
ordering given by Eq. (21). For example, for $k=3$ and $N=2$, $\mathcal{B}$ is
ordered as

$\mathcal{B}=\left(\left|0,0,2\right\rangle,\left|0,1,1\right\rangle,\left|0,2,0\right\rangle,\left|1,0,1\right\rangle,\left|1,1,0\right\rangle,\left|2,0,0\right\rangle\right).$
(22)

If the components of eigenstates of Bose-Hubbard chains would be independent,
their estimated covariance matrix would be close to the identity matrix.
Figure 9 shows that, even for $k=6$, there are significant off-diagonal terms,
indicating that the eigenstate components are not independent. The deviation
from identity matrix suggests that the mid-spectrum states deviate from
infinite-temperature states, even for significantly chaotic states, consistent
with Refs. Beugeling  _et al._ (2018); Luitz and Bar Lev (2016); Bäcker  _et
al._ (2019); Luitz  _et al._ (2020); Beugeling  _et al._ (2015b); Liu  _et
al._ (2018); LeBlond  _et al._ (2019); Haque  _et al._ (2020); Khaymovich  _et
al._ (2019). For small $k$, nonzero off-diagonal elements are even more
pronounced, and these appear in intricate patterns (Figure 9(a)).

##  Appendix D Variance and Covariance of EEVs of Random Eigenstates

In this section, we provide derivations of some of the equations that were
announced and used in the main text without proof. First, we derive the
results announced in Section III for random eigenstates with independent and
identically distributed (i.i.d.) components (Subsection D.1). Next, we prove
the result Eq. (20), announced in Subsection IV.3, for random vectors whose
components are not independent, with 2-point correlations described by a
covariance matrix (Subsection D.2).

###  D.1 Random vectors with i.i.d. components

Let $A$ be a $D\times D$ square matrix and $Z$ be a $D$-dimensional,
multivariate random state with identically and independently distributed
components $Z_{i}$, each with mean 0. The statistical variance of $\langle
A\rangle=\langle Z|A|Z\rangle$ is given by

$\displaystyle\operatorname{var}(\langle A\rangle)$
$\displaystyle=\sum_{i,j,i^{\prime},j^{\prime}=1}^{D}A_{i,j}A_{i^{\prime},j^{\prime}}\operatorname{cov}(Z_{i}Z_{j},Z_{i^{\prime}}Z_{j^{\prime}}),$
(23)

where

$\displaystyle\operatorname{cov}(Z_{i}Z_{j},Z_{i^{\prime}}Z_{j^{\prime}})$
$\displaystyle=E[Z_{i}Z_{j}Z_{i^{\prime}}Z_{j^{\prime}}]-E[Z_{i}Z_{j}]\cdot
E[Z_{i^{\prime}}Z_{j^{\prime}}]$
$\displaystyle=E[Z_{i}Z_{j}Z_{i^{\prime}}Z_{j^{\prime}}]-\delta_{ij}\delta_{i^{\prime}j^{\prime}}E[Z_{i}^{2}]^{2}.$
(24)

By the independence of $Z_{i},Z_{j}$ for $i\neq j$,
$E[Z_{i}Z_{j}Z_{i^{\prime}}Z_{j^{\prime}}]$ is only non-zero if there is no
index $i,j,i^{\prime},j^{\prime}$ different to the other three. The only
possibilities for this are

$\displaystyle i=j=i^{\prime}=j^{\prime}$ $\displaystyle i=j\text{ and
}i^{\prime}=j^{\prime}$ $\displaystyle i=i^{\prime}\text{ and }j=j^{\prime}$
$\displaystyle i=j^{\prime}\text{ and }j=i^{\prime},$ (25)

so

$\displaystyle E[Z_{i}$ $\displaystyle
Z_{j}Z_{i^{\prime}}Z_{j^{\prime}}]=\delta_{ij}\delta_{ii^{\prime}}\delta_{ij^{\prime}}E[Z_{i}^{4}]$
$\displaystyle+(1-\delta_{ij}\delta_{ii^{\prime}}\delta_{ij^{\prime}})$
$\displaystyle\times\left[\delta_{ij}\delta_{i^{\prime}j^{\prime}}E[Z_{i}^{2}]^{2}+\delta_{ii^{\prime}}\delta_{jj^{\prime}}E[Z_{i}^{2}]^{2}+\delta_{ij^{\prime}}\delta_{ji^{\prime}}E[Z_{i}^{2}]^{2}\right].$
(26)

Then we get

$\displaystyle\operatorname{var}(\langle A\rangle)$
$\displaystyle=E[Z_{1}^{4}]\sum_{i}A_{ii}^{2}+E[Z_{1}^{2}]^{2}$
$\displaystyle\times\left[\sum_{i\neq
i^{\prime}}A_{ii}A_{i^{\prime}i^{\prime}}+\sum_{i\neq j}A_{ij}^{2}+\sum_{i\neq
j}A_{ij}A_{ji}\right]$
$\displaystyle-E[Z_{1}^{2}]^{2}\sum_{i,i^{\prime}}A_{ii}A_{i^{\prime}i^{\prime}}$
(27)

and

$\displaystyle\operatorname{var}(\langle A\rangle)$
$\displaystyle=\left(E[Z_{1}^{4}]-E[Z_{1}^{2}]^{2}\right)\sum_{i}A_{ii}^{2}$
$\displaystyle+E[Z_{1}^{2}]^{2}\left[\sum_{i\neq j}A_{ij}^{2}+\sum_{i\neq
j}A_{ij}A_{ji}\right]$
$\displaystyle=\left(E[Z_{1}^{4}]-3E[Z_{1}^{2}]^{2}\right)\sum_{i}A_{ii}^{2}$
$\displaystyle+E[Z_{1}^{2}]^{2}\left[\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})\right].$
(28)

This is Eq. (III.1.1) of the main text. It is a rather general result, not
assuming a particular distribution of the components, only that they should be
independent and identically distributed and have mean 0.

Now, we specialize to the case that the components $Z_{i}$ are each normally
distributed with mean 0 and variance $1/D$. Then

$E[Z_{i}^{2}]=\frac{1}{D}E[(\sqrt{D}Z_{i})^{2}]=\frac{1}{D}$ (29)

and

$E[Z_{i}^{4}]=\frac{1}{D^{2}}E[(\sqrt{D}Z_{i})^{4}]=\frac{3}{D^{2}},$ (30)

because $\sqrt{DZ_{i}}$ is normally distributed with mean 0 and variance 1.
Plugging this into the above formula, we get

$\displaystyle\operatorname{var}(\langle A\rangle)$
$\displaystyle=\left(\frac{3}{D^{2}}-\frac{3}{D^{2}}\right)\sum_{i}A_{ii}^{2}$
(31)
$\displaystyle+\frac{1}{D^{2}}\left[\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})\right]$
$\displaystyle=\frac{1}{D^{2}}\left[\operatorname{tr}(A^{2})+\operatorname{tr}(AA^{\dagger})\right].$
(32)

This concludes the proof of Eq. (7).

###  D.2 Including 2-point correlations

We now relax the assumption of independence and allow the wavefunction
components to have 2-point correlations. The trace expression in this case is
Eq. (20), which we will now prove.

Assume $Z=LX$, where $X$ is a vector whose independent components are
Gaussian-distributed with mean zero and variance $1$, and $L$ is the Cholesky
root of a non-random $D\times D$ (covariance-) matrix $\Sigma$, i.e
$\Sigma=LL^{t}$. Then $Z$ has mean zero and the covariances between components
of $Z$ are given by $\Sigma$. We compute the statistical variance of the EEVs
as

$\displaystyle\operatorname{var}(\langle A\rangle)$
$\displaystyle=\operatorname{var}((LX)^{t}ALX)$
$\displaystyle=\operatorname{var}(X^{t}L^{t}ALX)$
$\displaystyle=\operatorname{var}(1/\sqrt{D}X^{t}DL^{t}AL1/\sqrt{D}X)$ (33)

Because the components of $1/\sqrt{D}X$ are normally distributed with mean 0
and variance $1/D$, we can use Eq. (7) with $A$ replaced by $L^{t}AL$. This
leads to

$\displaystyle\operatorname{var}(\langle A\rangle)$
$\displaystyle=\frac{1}{D^{2}}(\operatorname{tr}(D^{2}(L^{t}AL)^{2})+\operatorname{tr}(D^{2}L^{t}AL(L^{t}AL)^{t})$
$\displaystyle=\frac{1}{D^{2}}(\operatorname{tr}((D\Sigma
A)^{2})+\operatorname{tr}(D\Sigma AD\Sigma A^{t}),$ (34)

which concludes the proof of Eq. (20).

##  Appendix E Scaling of operators

In this section we derive results announced in the main text concerning the
scaling of the traces. For $A=a_{j}^{\dagger}a_{i}$ we prove the trace
expression of Eq. (8). For the classical limit, we prove the scaling relation,
Eq. (9), for $A$ being a (finite) linear combination of terms of the type
$a_{j}^{\dagger}a_{i}$.

Considering $A$ to be such a linear combination, $A^{2}$ and $A^{\dagger}A$
can be written as a sum of terms
$a_{j}^{\dagger}a_{i}a_{j^{\prime}}^{\dagger}a_{i^{\prime}}$. Consider a
defining basis state $|\psi\rangle=|n_{1},\dots,n_{k}\rangle$, which has
$n_{j}$ particles on site $j$. Then

$\displaystyle a_{j}^{\dagger}a_{i}|\psi\rangle$
$\displaystyle=\delta_{ij}n_{i}|\psi\rangle$ (35)
$\displaystyle+(1-\delta_{ij})\sqrt{n_{j}+1}\sqrt{n_{i}}|\dots,n_{i}-1,n_{j}+1,\dots\rangle.$

Using Eq. (35) twice we get

$\displaystyle\langle\psi|a_{j}^{\dagger}a_{i}$ $\displaystyle
a_{j^{\prime}}^{\dagger}a_{i^{\prime}}|\psi\rangle=\delta_{i^{\prime}j^{\prime}}\delta_{ij}n_{i^{\prime}}n_{i}$
$\displaystyle+(1-\delta_{i^{\prime}j^{\prime}})(1-\delta_{ij})\delta_{ij^{\prime}}\delta_{ji^{\prime}}(n_{j^{\prime}}+1)n_{i^{\prime}}.$
(36)

Now we will calculate
$\operatorname{tr}(a_{j}^{\dagger}a_{i}a_{j^{\prime}}^{\dagger}a_{i^{\prime}})$.
First let $i=j=i^{\prime}=j^{\prime}$. Using that there are
$\binom{N-l+k-2}{k-2}$ states with $l$ particles on site $i$, we get

$\operatorname{tr}(n_{i}^{2})=\sum_{l=0}^{N}l^{2}\binom{N-l+k-2}{k-2}.$ (37)

Writing $l^{2}$ in terms of binomial coefficients and exploiting an upper
index Vandermonde identity, namely

$\sum_{l=0}^{n}\binom{l}{c_{1}}\binom{n-l}{c_{2}}=\binom{n+1}{c_{1}+c_{2}+1},$
(38)

for constants $n$, $c_{1}$ and $c_{2}$, we get

$\operatorname{tr}(n_{i}^{2})=\frac{N(2N+k-1)}{k(k+1)}D.$ (39)

Now let $i=j$ and $i^{\prime}=j^{\prime}$ but $i\neq i^{\prime}$. There are
$\binom{N-l-s+k-3}{k-3}$ many states with $l$ particles on site $i$ and $s$
particles on site $i^{\prime}$, so

$\operatorname{tr}(n_{i}n_{i^{\prime}})=\sum_{l=0}^{N}\sum_{s=0}^{N-l}ls\binom{N-l-s+k-3}{k-3}.$
(40)

Invoking the Vandermonde identity Eq. (38) twice gives us

$\operatorname{tr}(n_{i}n_{i^{\prime}})=\frac{N(N-1)}{(k+1)k}D.$ (41)

The case $i\neq j$ and $i^{\prime}\neq j^{\prime}$, but $j=i^{\prime}$ and
$i=j^{\prime}$, works exactly the same as $i=j$ and $i^{\prime}=j^{\prime}$
but $i\neq i^{\prime}$. Using Eq. (38) twice on

$\operatorname{tr}(a_{j}^{\dagger}a_{i}a_{i}^{\dagger}a_{j})=\sum_{l=0}^{N}\sum_{s=0}^{N-l}(l+1)s\binom{N-l-s+k-3}{k-3}$
(42)

results in

$\operatorname{tr}(a_{j}^{\dagger}a_{i}a_{i}^{\dagger}a_{j})=\frac{N(N+k)}{k(k+1)}D.$
(43)

In the classical limit, $k\ll N$, Eqs. (39), (41) and (43) all scale as $\sim
N^{2}D$, so the corresponding trace expressions of all typical observables
scale as

$\operatorname{tr}(A^{2})+\operatorname{tr}(A^{\dagger}A)\sim N^{2}D.$ (44)

In the thermodynamic limit, $k=\rho N$ for a density $\rho$ not depending on
$N$ nor $k$ and $k,N\to\infty$, the trace expressions scale as

$\operatorname{tr}(A^{2})+\operatorname{tr}(A^{\dagger}A)\sim D,$ (45)

because the leading order of $k$ and $N$ is quadratic in both the numerators
and the denumerators in Eqs. (39), (41) and (43). By combining this with the
results in Eq. (D) we rediscover the $D^{-1/2}$ dependence of the EEV
fluctuations in the thermodynamic limit for Gaussian states and observables
$A$, where $\operatorname{tr}(A^{2})$ and $\operatorname{tr}(AA^{\dagger})$
are bounded in $D$.

## References

  * Deutsch (1991) J. M. Deutsch, “Quantum statistical mechanics in a closed system,” Phys. Rev. A 43, 2046–2049 (1991).
  * Srednicki (1994) Mark Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
  * Srednicki (1996) Mark Srednicki, “Thermal fluctuations in quantized chaotic systems,” Journal of Physics A: Mathematical and General 29, L75 (1996).
  * Srednicki (1999) Mark Srednicki, “The approach to thermal equilibrium in quantized chaotic systems,” Journal of Physics A: Mathematical and General 32, 1163 (1999).
  * Rigol  _et al._ (2008) Marcos Rigol, Vanja Dunjko, and Maxim Olshanii, “Thermalization and its mechanism for generic isolated quantum systems,” Nature 452, 854–858 (2008).
  * Polkovnikov  _et al._ (2011) Anatoli Polkovnikov, Krishnendu Sengupta, Alessandro Silva, and Mukund Vengalattore, “Colloquium: Nonequilibrium dynamics of closed interacting quantum systems,” Rev. Mod. Phys. 83, 863–883 (2011).
  * Rigol and Srednicki (2012) Marcos Rigol and Mark Srednicki, “Alternatives to eigenstate thermalization,” Phys. Rev. Lett. 108, 110601 (2012).
  * Eisert  _et al._ (2015) Jens Eisert, Mathis Friesdorf, and Christian Gogolin, “Quantum many-body systems out of equilibrium,” Nature Physics 11, 124–130 (2015).
  * D’Alessio  _et al._ (2016) Luca D’Alessio, Yariv Kafri, Anatoli Polkovnikov, and Marcos Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65, 239 (2016).
  * Neuenhahn and Marquardt (2012) Clemens Neuenhahn and Florian Marquardt, “Thermalization of interacting fermions and delocalization in fock space,” Phys. Rev. E 85, 060101 (2012).
  * Beugeling  _et al._ (2014) W. Beugeling, R. Moessner, and Masudul Haque, “Finite-size scaling of eigenstate thermalization,” Phys. Rev. E 89, 042112 (2014).
  * Beugeling  _et al._ (2015a) Wouter Beugeling, Roderich Moessner, and Masudul Haque, “Off-diagonal matrix elements of local operators in many-body quantum systems,” Phys. Rev. E 91, 012144 (2015a).
  * Ziraldo and Santoro (2013) Simone Ziraldo and Giuseppe E. Santoro, “Relaxation and thermalization after a quantum quench: Why localization is important,” Phys. Rev. B 87, 064201 (2013).
  * Ikeda  _et al._ (2013) Tatsuhiko N. Ikeda, Yu Watanabe, and Masahito Ueda, “Finite-size scaling analysis of the eigenstate thermalization hypothesis in a one-dimensional interacting bose gas,” Phys. Rev. E 87, 012125 (2013).
  * Alba (2015) Vincenzo Alba, “Eigenstate thermalization hypothesis and integrability in quantum spin chains,” Phys. Rev. B 91, 155123 (2015).
  * Nandy  _et al._ (2016) Sourav Nandy, Arnab Sen, Arnab Das, and Abhishek Dhar, “Eigenstate gibbs ensemble in integrable quantum systems,” Phys. Rev. B 94, 245131 (2016).
  * Magán (2016) Javier M. Magán, “Random free fermions: An analytical example of eigenstate thermalization,” Phys. Rev. Lett. 116, 030401 (2016).
  * Haque and McClarty (2019) M. Haque and P. A. McClarty, “Eigenstate thermalization scaling in majorana clusters: From chaotic to integrable sachdev-ye-kitaev models,” Phys. Rev. B 100, 115122 (2019).
  * Mierzejewski and Vidmar (2020) Marcin Mierzejewski and Lev Vidmar, “Quantitative impact of integrals of motion on the eigenstate thermalization hypothesis,” Phys. Rev. Lett. 124, 040603 (2020).
  * LeBlond  _et al._ (2019) Tyler LeBlond, Krishnanand Mallayya, Lev Vidmar, and Marcos Rigol, “Entanglement and matrix elements of observables in interacting integrable systems,” Phys. Rev. E 100, 062134 (2019).
  * Rigol (2009) Marcos Rigol, “Breakdown of thermalization in finite one-dimensional systems,” Phys. Rev. Lett. 103, 100403 (2009).
  * Biroli  _et al._ (2010) Giulio Biroli, Corinna Kollath, and Andreas M. Läuchli, “Effect of rare fluctuations on the thermalization of isolated quantum systems,” Phys. Rev. Lett. 105, 250401 (2010).
  * Rigol and Santos (2010) Marcos Rigol and Lea F. Santos, “Quantum chaos and thermalization in gapped systems,” Phys. Rev. A 82, 011604 (2010).
  * Santos and Rigol (2010) Lea F. Santos and Marcos Rigol, “Localization and the effects of symmetries in the thermalization properties of one-dimensional quantum systems,” Phys. Rev. E 82, 031130 (2010).
  * Roux (2010) Guillaume Roux, “Finite-size effects in global quantum quenches: Examples from free bosons in an harmonic trap and the one-dimensional Bose-Hubbard model,” Phys. Rev. A 81, 053604 (2010).
  * Motohashi (2011) Atsushi Motohashi, “Thermalization of atom-molecule bose gases in a double-well potential,” Phys. Rev. A 84, 063631 (2011).
  * Brandino  _et al._ (2012) G. P. Brandino, A. De Luca, R. M. Konik, and G. Mussardo, “Quench dynamics in randomly generated extended quantum models,” Phys. Rev. B 85, 214435 (2012).
  * Steinigeweg  _et al._ (2013) R. Steinigeweg, J. Herbrych, and P. Prelovšek, “Eigenstate thermalization within isolated spin-chain systems,” Phys. Rev. E 87, 012118 (2013).
  * Kim  _et al._ (2014) Hyungwon Kim, Tatsuhiko N. Ikeda, and David A. Huse, “Testing whether all eigenstates obey the eigenstate thermalization hypothesis,” Phys. Rev. E 90, 052105 (2014).
  * Sorg  _et al._ (2014) S. Sorg, L. Vidmar, L. Pollet, and F. Heidrich-Meisner, “Relaxation and thermalization in the one-dimensional bose-hubbard model: A case study for the interaction quantum quench from the atomic limit,” Phys. Rev. A 90, 033606 (2014).
  * Steinigeweg  _et al._ (2014) R. Steinigeweg, A. Khodja, H. Niemeyer, C. Gogolin, and J. Gemmer, “Pushing the limits of the eigenstate thermalization hypothesis towards mesoscopic quantum systems,” Phys. Rev. Lett. 112, 130403 (2014).
  * Fratus and Srednicki (2015) Keith R. Fratus and Mark Srednicki, “Eigenstate thermalization in systems with spontaneously broken symmetry,” Phys. Rev. E 92, 040103 (2015).
  * Nandkishore and Huse (2015) Rahul Nandkishore and David A. Huse, “Many-Body Localization and Thermalization in Quantum Statistical Mechanics,” Annual Review of Condensed Matter Physics 6, 15–38 (2015).
  * Khodja  _et al._ (2015) Abdellah Khodja, Robin Steinigeweg, and Jochen Gemmer, “Relevance of the eigenstate thermalization hypothesis for thermal relaxation,” Phys. Rev. E 91, 012120 (2015).
  * Mondaini  _et al._ (2016) Rubem Mondaini, Keith R. Fratus, Mark Srednicki, and Marcos Rigol, “Eigenstate thermalization in the two-dimensional transverse field ising model,” Phys. Rev. E 93, 032104 (2016).
  * Chandran  _et al._ (2016) A. Chandran, Marc D. Schulz, and F. J. Burnell, “The eigenstate thermalization hypothesis in constrained hilbert spaces: A case study in non-abelian anyo n chains,” Phys. Rev. B 94, 235122 (2016).
  * Luitz and Bar Lev (2016) David J. Luitz and Yevgeny Bar Lev, “Anomalous thermalization in ergodic systems,” Phys. Rev. Lett. 117, 170404 (2016).
  * Mondaini and Rigol (2017) Rubem Mondaini and Marcos Rigol, “Eigenstate thermalization in the two-dimensional transverse field ising model. ii. off-diagonal matrix elements of observables,” Phys. Rev. E 96, 012157 (2017).
  * Sonner and Vielma (2017) Julian Sonner and Manuel Vielma, “Eigenstate thermalization in the sachdev-ye-kitaev model,” Journal of High Energy Physics 2017, 149 (2017).
  * Lan and Powell (2017) Zhihao Lan and Stephen Powell, “Eigenstate thermalization hypothesis in quantum dimer models,” Phys. Rev. B 96, 115140 (2017).
  * Yoshizawa  _et al._ (2018) Toru Yoshizawa, Eiki Iyoda, and Takahiro Sagawa, “Numerical large deviation analysis of the eigenstate thermalization hypothesis,” Phys. Rev. Lett. 120, 200604 (2018).
  * Dymarsky  _et al._ (2018) Anatoly Dymarsky, Nima Lashkari, and Hong Liu, “Subsystem eigenstate thermalization hypothesis,” Phys. Rev. E 97, 012140 (2018).
  * Hunter-Jones  _et al._ (2018) Nicholas Hunter-Jones, Junyu Liu, and Yehao Zhou, “On thermalization in the syk and supersymmetric syk models,” Journal of High Energy Physics 2018, 142 (2018).
  * Jansen  _et al._ (2019) David Jansen, Jan Stolpp, Lev Vidmar, and Fabian Heidrich-Meisner, “Eigenstate thermalization and quantum chaos in the holstein polaron model,” Phys. Rev. B 99, 155130 (2019).
  * Khaymovich  _et al._ (2019) Ivan M. Khaymovich, Masudul Haque, and Paul A. McClarty, “Eigenstate thermalization, random matrix theory, and behemoths,” Phys. Rev. Lett. 122, 070601 (2019).
  * Hamazaki and Ueda (2018) Ryusuke Hamazaki and Masahito Ueda, “Atypicality of most few-body observables,” Phys. Rev. Lett. 120, 080603 (2018).
  * Khemani  _et al._ (2019) Vedika Khemani, Chris R. Laumann, and Anushya Chandran, “Signatures of integrability in the dynamics of rydberg-blockaded chains,” Phys. Rev. B 99, 161101 (2019).
  * Brenes  _et al._ (2020) Marlon Brenes, Tyler LeBlond, John Goold, and Marcos Rigol, “Eigenstate thermalization in a locally perturbed integrable system,” Phys. Rev. Lett. 125, 070605 (2020).
  * Polkovnikov (2003) Anatoli Polkovnikov, “Quantum corrections to the dynamics of interacting bosons: Beyond the truncated wigner approximation,” Phys. Rev. A 68, 053604 (2003).
  * Mahmud  _et al._ (2005) Khan W. Mahmud, Heidi Perry, and William P. Reinhardt, “Quantum phase-space picture of bose-einstein condensates in a double well,” Phys. Rev. A 71, 023615 (2005).
  * Hiller  _et al._ (2006) Moritz Hiller, Tsampikos Kottos, and T. Geisel, “Complexity in parametric bose-hubbard hamiltonians and structural analysis of eigenstates,” Phys. Rev. A 73, 061604 (2006).
  * Mossmann and Jung (2006) S. Mossmann and C. Jung, “Semiclassical approach to bose-einstein condensates in a triple well potential,” Phys. Rev. A 74, 033601 (2006).
  * Graefe and Korsch (2007) E. M. Graefe and H. J. Korsch, “Semiclassical quantization of an $n$-particle bose-hubbard model,” Phys. Rev. A 76, 032116 (2007).
  * Trimborn  _et al._ (2009) F. Trimborn, D. Witthaut, and H. J. Korsch, “Beyond mean-field dynamics of small bose-hubbard systems based on the number-conserving phase-space approach,” Phys. Rev. A 79, 013608 (2009).
  * Polkovnikov (2010) Anatoli Polkovnikov, “Phase space representation of quantum dynamics,” Annals of Physics 325, 1790 (2010).
  * Pawłowski  _et al._ (2011) Krzysztof Pawłowski, Paweł Ziń, Kazimierz Rzażewski, and Marek Trippenbach, “Revivals in an attractive bose-einstein condensate in a double-well potential and their decoherence,” Phys. Rev. A 83, 033606 (2011).
  * Simon and Strunz (2012) Lena Simon and Walter T. Strunz, “Analytical results for josephson dynamics of ultracold bosons,” Phys. Rev. A 86, 053625 (2012).
  * Simon and Strunz (2014) Lena Simon and Walter T. Strunz, “Time-dependent semiclassics for ultracold bosons,” Phys. Rev. A 89, 052112 (2014).
  * Veksler and Fishman (2015) Hagar Veksler and Shmuel Fishman, “Semiclassical analysis of bose–hubbard dynamics,” New Journal of Physics 17, 053030 (2015).
  * Engl  _et al._ (2015) Thomas Engl, Juan Diego Urbina, and Klaus Richter, “Periodic mean-field solutions and the spectra of discrete bosonic fields: Trace formula for bose-hubbard models,” Phys. Rev. E 92, 062907 (2015).
  * Engl  _et al._ (2016) Thomas Engl, Juan Diego Urbina, and Klaus Richter, “The semiclassical propagator in fock space: dynamical echo and many-body interference,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150159 (2016).
  * Ray  _et al._ (2016a) Shouryya Ray, Paula Ostmann, Lena Simon, Frank Grossmann, and Walter T. Strunz, “Dynamics of interacting bosons using the herman-kluk semiclassical initial value representation,” Journal of Physics A: Mathematical and Theoretical 49, 165303 (2016a).
  * Dubertrand and Müller (2016) Rémy Dubertrand and Sebastian Müller, “Spectral statistics of chaotic many-body systems,” New Journal of Physics 18, 033009 (2016).
  * Tomsovic  _et al._ (2018) Steven Tomsovic, Peter Schlagheck, Denis Ullmo, Juan-Diego Urbina, and Klaus Richter, “Post-ehrenfest many-body quantum interferences in ultracold atoms far out of equilibrium,” Phys. Rev. A 97, 061606 (2018).
  * Tomsovic (2018) Steven Tomsovic, “Complex saddle trajectories for multidimensional quantum wave packet and coherent state propagation: Application to a many-body system,” Phys. Rev. E 98, 023301 (2018).
  * Rammensee  _et al._ (2018) Josef Rammensee, Juan Diego Urbina, and Klaus Richter, “Many-Body Quantum Interference and the Saturation of Out-of-Time-Order Correlators,” Physical Review Letters 121, 124101 (2018).
  * Kidd  _et al._ (2019) R. A. Kidd, M. K. Olsen, and J. F. Corney, “Quantum chaos in a bose-hubbard dimer with modulated tunneling,” Phys. Rev. A 100, 013625 (2019).
  * Schlagheck  _et al._ (2019) Peter Schlagheck, Denis Ullmo, Juan Diego Urbina, Klaus Richter, and Steven Tomsovic, “Enhancement of many-body quantum interference in chaotic bosonic systems: The role of symmetry and dynamics,” Phys. Rev. Lett. 123, 215302 (2019).
  * Nemoto  _et al._ (2000) K. Nemoto, C. A. Holmes, G. J. Milburn, and W. J. Munro, “Quantum dynamics of three coupled atomic bose-einstein condensates,” Phys. Rev. A 63, 013604 (2000).
  * Weiss and Teichmann (2008) Christoph Weiss and Niklas Teichmann, “Differences between mean-field dynamics and $n$-particle quantum dynamics as a signature of entanglement,” Phys. Rev. Lett. 100, 140408 (2008).
  * Viscondi and Furuya (2011) Thiago F Viscondi and K Furuya, “Dynamics of a bose–einstein condensate in a symmetric triple-well trap,” Journal of Physics A: Mathematical and Theoretical 44, 175301 (2011).
  * Gertjerenken and Weiss (2013) Bettina Gertjerenken and Christoph Weiss, “Beyond-mean-field behavior of large bose-einstein condensates in double-well potentials,” Phys. Rev. A 88, 033608 (2013).
  * Kolovsky (2016) Andrey R. Kolovsky, “Bose-hubbard hamiltonian: Quantum chaos approach,” Int. J. Mod. Phys. B 30, 1630009 (2016).
  * Heinisch and Holthaus (2016) Christoph Heinisch and Martin Holthaus, “Entropy production within a pulsed bose–einstein condensate,” Zeitschrift für Naturforschung A 71, 875 (2016).
  * Rautenberg and Gärttner (2020) Michael Rautenberg and Martin Gärttner, “Classical and quantum chaos in a three-mode bosonic system,” Phys. Rev. A 101, 053604 (2020).
  * Hiller  _et al._ (2009) Moritz Hiller, Tsampikos Kottos, and Theo Geisel, “Wave-packet dynamics in energy space of a chaotic trimeric bose-hubbard system,” Phys. Rev. A 79, 023621 (2009).
  * Arwas  _et al._ (2014) Geva Arwas, Amichay Vardi, and Doron Cohen, “Triangular bose-hubbard trimer as a minimal model for a superfluid circuit,” Phys. Rev. A 89, 013601 (2014).
  * Han and Wu (2016) Xizhi Han and Biao Wu, “Ehrenfest breakdown of the mean-field dynamics of bose gases,” Phys. Rev. A 93, 023621 (2016).
  * Bürkle and Anglin (2019) R. Bürkle and J. R. Anglin, “Threshold coupling strength for equilibration between small systems,” Phys. Rev. A 99, 063617 (2019).
  * Bychek  _et al._ (2020) Anna A. Bychek, Pavel S. Muraev, Dmitrii N. Maksimov, Evgeny N. Bulgakov, and Andrey R. Kolovsky, “Chaotic and regular dynamics in the three-site bose-hubbard model,” AIP Conference Proceedings 2241, 020007 (2020).
  * Kolovsky and Buchleitner (2004) A. R Kolovsky and A Buchleitner, “Quantum chaos in the bose-hubbard model,” Europhysics Letters (EPL) 68, 632 (2004).
  * Kollath  _et al._ (2010) Corinna Kollath, Guillaume Roux, Giulio Biroli, and Andreas M Läuchli, “Statistical properties of the spectrum of the extended bose–hubbard model,” Journal of Statistical Mechanics: Theory and Experiment 2010, P08011 (2010).
  * de la Cruz  _et al._ (2020) Javier de la Cruz, Sergio Lerma-Hernández, and Jorge G. Hirsch, “Quantum chaos in a system with high degree of symmetries,” Phys. Rev. E 102, 032208 (2020).
  * Pausch  _et al._ (2020) Lukas Pausch, Edoardo G Carnio, Alberto Rodríguez, and Andreas Buchleitner, “Chaos and ergodicity across the energy spectrum of interacting bosons,” arXiv preprint arXiv:2009.05295 (2020).
  * Beugeling  _et al._ (2018) Wouter Beugeling, Arnd Bäcker, Roderich Moessner, and Masudul Haque, “Statistical properties of eigenstate amplitudes in complex quantum systems,” Phys. Rev. E 98, 022204 (2018).
  * Bäcker  _et al._ (2019) Arnd Bäcker, Masudul Haque, and Ivan M Khaymovich, “Multifractal dimensions for random matrices, chaotic quantum maps, and many-body systems,” Phys. Rev. E 100, 032117 (2019).
  * Luitz  _et al._ (2020) David J. Luitz, Ivan M. Khaymovich, and Yevgeny Bar Lev, “Multifractality and its role in anomalous transport in the disordered XXZ spin-chain,” SciPost Phys. Core 2, 6 (2020).
  * Beugeling  _et al._ (2015b) W Beugeling, A Andreanov, and Masudul Haque, “Global characteristics of all eigenstates of local many-body hamiltonians: participation ratio and entanglement entropy,” Journal of Statistical Mechanics: Theory and Experiment 2015, P02002 (2015b).
  * Liu  _et al._ (2018) Chunxiao Liu, Xiao Chen, and Leon Balents, “Quantum entanglement of the sachdev-ye-kitaev models,” Phys. Rev. B 97, 245126 (2018).
  * Haque  _et al._ (2020) Masudul Haque, Paul A McClarty, and Ivan M Khaymovich, “Entanglement of mid-spectrum eigenstates of chaotic many-body systems–deviation from random ensembles,” arXiv preprint arXiv:2008.12782 (2020).
  * Gemmer  _et al._ (2009) J. Gemmer, M. Michel, and G. Mahler,  _Quantum Thermodynamics: Emergence of Thermodynamic Behavior Within Composite Quantum Systems_, Lecture Notes in Physics (Springer, 2009).
  * Reimann (2007) Peter Reimann, “Typicality for generalized microcanonical ensembles,” Phys. Rev. Lett. 99, 160404 (2007).
  * Reimann (2008) Peter Reimann, “Typicality of pure states randomly sampled according to the gaussian adjusted projected measure,” Journal of Statistical Physics 132, 921 (2008).
  * Lloyd (2013) Seth Lloyd, “Pure state quantum statistical mechanics and black holes,” arXiv preprint arXiv:1307.0378 (2013).
  * Feingold and Peres (1983) Mario Feingold and Asher Peres, “Regular and chaotic motion of coupled rotators,” Physica D: Nonlinear Phenomena 9, 433 (1983).
  * Feingold  _et al._ (1984) Mario Feingold, Nimrod Moiseyev, and Asher Peres, “Ergodicity and mixing in quantum theory. ii,” Phys. Rev. A 30, 509 (1984).
  * Robb and Reichl (1998) Daniel T. Robb and L. E. Reichl, “Chaos in a two-spin system with applied magnetic field,” Phys. Rev. E 57, 2458–2459 (1998).
  * Emerson and Ballentine (2001) J. Emerson and L.E. Ballentine, “Characteristics of quantum-classical correspondence for two interacting spins,” Phys. Rev. A 63, 052103 (2001).
  * Ballentine (2004) L. E. Ballentine, “Quantum-to-classical limit in a hamiltonian system,” Phys. Rev. A 70, 032111 (2004).
  * Fan  _et al._ (2017) Yiyun Fan, Sven Gnutzmann, and Yuqi Liang, “Quantum chaos for nonstandard symmetry classes in the feingold-peres model of coupled tops,” Phys. Rev. E 96, 062207 (2017).
  * Pappalardi  _et al._ (2018) Silvia Pappalardi, Angelo Russomanno, Bojan Žunkovič, Fernando Iemini, Alessandro Silva, and Rosario Fazio, “Scrambling and entanglement spreading in long-range spin chains,” Phys. Rev. B 98, 134303 (2018).
  * Piga  _et al._ (2019) Angelo Piga, Maciej Lewenstein, and James Q. Quach, “Quantum chaos and entanglement in ergodic and nonergodic systems,” Phys. Rev. E 99, 032213 (2019).
  * Lerose and Pappalardi (2020) Alessio Lerose and Silvia Pappalardi, “Bridging entanglement dynamics and chaos in semiclassical systems,” Phys. Rev. A 102, 032404 (2020).
  * Xu  _et al._ (2020) Tianrui Xu, Thomas Scaffidi, and Xiangyu Cao, “Does scrambling equal chaos?” Phys. Rev. Lett. 124, 140602 (2020).
  * Pappalardi  _et al._ (2020) Silvia Pappalardi, Anatoli Polkovnikov, and Alessandro Silva, “Quantum echo dynamics in the Sherrington-Kirkpatrick model,” SciPost Phys. 9, 21 (2020).
  * Mondal  _et al._ (2020) Debabrata Mondal, Sudip Sinha, and S. Sinha, “Chaos and quantum scars in a coupled top model,” Phys. Rev. E 102, 020101 (2020).
  * Bastarrachea-Magnani  _et al._ (2014) M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch, “Comparative quantum and semiclassical analysis of atom-field systems. ii. chaos and regularity,” Phys. Rev. A 89, 032102 (2014).
  * Bastarrachea-Magnani  _et al._ (2015) Miguel Angel Bastarrachea-Magnani, Baldemar López del Carpio, Sergio Lerma-Hernández, and Jorge G Hirsch, “Chaos in the dicke model: quantum and semiclassical analysis,” Physica Scripta 90, 068015 (2015).
  * Ray  _et al._ (2016b) S. Ray, A. Ghosh, and S. Sinha, “Quantum signature of chaos and thermalization in the kicked dicke model,” Phys. Rev. E 94, 032103 (2016b).
  * Bastarrachea-Magnani  _et al._ (2016) M. A. Bastarrachea-Magnani, B. López-del Carpio, J. Chávez-Carlos, S. Lerma-Hernández, and J. G. Hirsch, “Delocalization and quantum chaos in atom-field systems,” Phys. Rev. E 93, 022215 (2016).
  * Lewis-Swan  _et al._ (2019) RJ Lewis-Swan, Arghavan Safavi-Naini, John J Bollinger, and Ana M Rey, “Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the dicke model,” Nature communications 10, 1–9 (2019).
  * Chávez-Carlos  _et al._ (2019) Jorge Chávez-Carlos, B. López-del Carpio, Miguel A. Bastarrachea-Magnani, Pavel Stránský, Sergio Lerma-Hernández, Lea F. Santos, and Jorge G. Hirsch, “Quantum and classical lyapunov exponents in atom-field interaction systems,” Phys. Rev. Lett. 122, 024101 (2019).
  * Wang and Robnik (2020) Qian Wang and Marko Robnik, “Statistical properties of the localization measure of chaotic eigenstates in the dicke model,” Phys. Rev. E 102, 032212 (2020).
  * Villaseñor  _et al._ (2020) D Villaseñor, S Pilatowsky-Cameo, M A Bastarrachea-Magnani, S Lerma-Hernández, L F Santos, and J G Hirsch, “Quantum vs classical dynamics in a spin-boson system: manifestations of spectral correlations and scarring,” New Journal of Physics 22, 063036 (2020).
  * Pilatowsky-Cameo  _et al._ (2020) Saúl Pilatowsky-Cameo, Jorge Chávez-Carlos, Miguel A. Bastarrachea-Magnani, Pavel Stránský, Sergio Lerma-Hernández, Lea F. Santos, and Jorge G. Hirsch, “Positive quantum lyapunov exponents in experimental systems with a regular classical limit,” Phys. Rev. E 101, 010202 (2020).
  * Sinha and Sinha (2020) Sudip Sinha and S. Sinha, “Chaos and quantum scars in bose-josephson junction coupled to a bosonic mode,” Phys. Rev. Lett. 125, 134101 (2020).
  * Oganesyan and Huse (2007) Vadim Oganesyan and David A. Huse, “Localization of interacting fermions at high temperature,” Phys. Rev. B 75, 155111 (2007).
  * Atas  _et al._ (2013) Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, “Distribution of the ratio of consecutive level spacings in random matrix ensembles,” Phys. Rev. Lett. 110, 084101 (2013).