Ergodicity probes: using time-fluctuations to measure the
Hilbert space dimension
Charlie Nation
1
and Diego Porras
2
1
Department of Physics and Astronomy, University of Sussex, Brighton, BN1 9QH, United Kingdom.
2
Institute of Fundamental Physics, IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain
Quantum devices, such as quantum sim-
ulators, quantum annealers, and quan-
tum computers, may be exploited to solve
problems beyond what is tractable with
classical computers. This may be achieved
as the Hilbert space available to perform
such ‘calculations’ is far larger than that
which may be classically simulated. In
practice, however, quantum devices have
imperfections, which may limit the accessi-
bility to the whole Hilbert space. We thus
determine that the dimension of the space
of quantum states that are available to a
quantum device is a meaningful measure of
its functionality, though unfortunately this
quantity cannot be directly experimentally
determined. Here we outline an experi-
mentally realisable approach to obtaining
the required Hilbert space dimension of
such a device to compute its time evolu-
tion, by exploiting the thermalization dy-
namics of a probe qubit. This is achieved
by obtaining a fluctuation-dissipation the-
orem for high-temperature chaotic quan-
tum systems, which facilitates the extrac-
tion of information on the Hilbert space
dimension via measurements of the decay
rate, and time-fluctuations.
1 Introduction
The ability to control and manipulate micro-
scopic systems at the single particle level is an es-
sential requirement for many quantum technolo-
gies. Experimental setups where atoms or qubits
can be arranged in ordered structures and studied
in quantum non-equilibrium states include neu-
tral atoms in optical lattices [1–3], trapped ions
Charlie Nation: C.Nation@sussex.ac.uk
Diego Porras: D.Porras@iff.csic.es
[4–7], Rydberg atoms [8, 9], and superconducting
circuits[10, 11]. These systems can be used for
the quantum simulation of many-body models,
or different forms of digital or adiabatic quantum
computing. Most of these physical setups have
limitations in the accessibility to certain observ-
ables. Thus, having extra tools to characterize
quantum systems in a simple and efficient way
can be useful in the diagnosis and certification of
quantum devices.
One of the most prominent properties of a
quantum device is its size in terms of the di-
mension of the associated Hilbert space. The
size of a quantum computer or simulator is of-
ten given in terms of number of qubits, such that
the Hilbert space dimension is 2
N
for N qubits.
This is a measure that ignores the effect of disor-
der or the possible lack of connectivity between
different zones in the device.
A more useful quantity would be the number
of eigenstates of the Hamiltonian that take part
in the quantum dynamics, which is bounded in
this case by 2
N
, but would exclude the degrees of
freedom that do not contribute to the evolution of
the initial state. Thereby establishing an effective
Hilbert space dimension that more accurately de-
scribes the complexity of the system, in terms of
some effective fully connected Hamiltonian. This
Hilbert space dimension is, however, an elusive
measure in realistic experimental situations.
In this work we show that the equilibration dy-
namics [12–18] of a quantum system can be used
to extract such information on the dimension of
the Hilbert space of a quantum device, in terms
of the effective number of states that contribute
to the dynamics of a local observable. Indeed,
advancements in quantum technologies described
above have inspired a bounty of theoretical work
in the field of quantum thermalization [19–30]. In
the following, we aim to help ‘bridge the gap’ be-
tween theoretical and experimental work in this
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 1
arXiv:1906.06206v3 [quant-ph] 19 Nov 2019
field [31].
We assume a quantum quench scenario [8, 21,
28, 32] in which a quantum system is initialized
in a fully-decohered, infinite temperature state,
except for a subsystem that acts as a sensor and
is prepared in a pure state. For simplicity, we as-
sume that this subsystem is a single qubit, which
we refer to as the ‘probe’ qubit. The relaxation
dynamics of the probe qubit depends on the de-
tails of the underlying structure of the Hamilto-
nian, however, in the most generic case of non-
integrable systems, an appropriate description
can be given in terms of random matrix theory
(RMT) [25, 30, 33–38]. We show that the time-
fluctuations of the probe in the long-time limit
contain information about the Hilbert space di-
mension of the device.
Our article is structured as follows. Firstly,
in Section 2 we present the basic scheme and
summarize our main result, which relies on an
infinite-temperature fluctuation dissipation the-
orem (FDT) [39] for the dynamics of the probe
qubit in order to extract information on the scal-
ing of the effective Hilbert space of a quantum de-
vice. We continue, presenting numerical calcula-
tions that validate our predictions via exact diag-
onalization of a spin chain Hamiltonian in Section
3. We then outline in more detail our RMT model
in Section 4. In Section 5 we derive an expression
for the time dependence of generic observables in
chaotic quantum systems, and discuss how this
can be exploited for our measurements. We then
derive the FDT, and extensions from the RMT
model, and to finite temperatures, in Section 6.
Finally, we summarise our findings in Section 7.
Further details and proofs are included in Appen-
dices. We note that this is arranged such that
our key findings can be understood from sections
2 and 3, with the detailed calculations presented
later in the text.
2 Setup and main results
2.1 Proposed setup
We assume that we have a quantum system (from
here on, ‘quantum device’) made up of a single
‘probe’ qubit, initialized in a pure state, which at
t = 0 is coupled to the rest of the quantum device
(referred to as the ‘bath’). We initialize the bath
in an infinite temperature state. This is sketched
in Figure 1. This can be routinely achieved, for
Figure 1: Illustration of our proposed setup (quantum
device). A single qubit, labelled ‘Probe’, is coupled lo-
cally to a part of a larger non-integrable bath, initialized
in an infinite temperature state ρ
B
∼
ˆ
I (this restriction
is removed below). Experimentally, for our protocol, one
needs access only to an observable of the Probe qubit.
example, in a quantum computer device that has
not been properly initialized. Since quantum de-
vices always suffer from some kind of decoher-
ence, creating an infinite temperature (or fully
decohered) state is typically a simple task. We
then let this qubit evolve in time and reach an
equilibrium state. The initial state is thus,
ρ(t = 0) = | ↑ih↑| ⊗ ρ
B
. (1)
ρ
B
is the density matrix of the bath in the fully
decohered state,
ρ
B
=
1
N
B
N
B
X
β
|φ
B,β
ihφ
B,β
|, (2)
where |φ
B,α
B
i is a set of orthonormal wave func-
tions in the bath, which we take as the eigenstates
of the bath Hamiltonian H
B
. N
B
is the Hilbert
space dimension of the bath, and is the param-
eter on which we wish to infer information via
measurements of the probe qubit. Later on the
fully decohered state condition will be relaxed,
allowing for finite temperatures.
The system evolves under the interacting
Hamiltonian,
H = H
0
+ V, (3)
where H
0
=
ˆ
1 ⊗ H
B
is the Hamiltonian of the
uncoupled probe + bath device, and we assume
that the probe qubit does not evolve at all in the
absence of coupling to the bath, which is given by
the operator V . A prominent role in the following
discussions will be played by the eigensystems of
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 2
the interacting Hamiltonian,
H|ψ
µ
i = E
µ
|ψ
µ
i, (4)
and the eigensystem of the uncoupled probe-bath
system
H
0
|φ
α
i = E
(0)
α
|φ
α
i. (5)
All throughout this work we will use the conven-
tion that indices µ, ν refer to summations over
eigenstates of H, whereas α, β refer to eigenstates
of H
0
.
We focus on the quantum dynamics of the ex-
pectation value of a probe observable, which we
take for concreteness to be hσ
z
(t)i = Tr(ρ(t)σ
z
).
Note that, however, the results should not de-
pend on the particular choice for the probe ob-
servable, since H
0
has no probe energy term. The
probe-bath coupling V , may introduce a depen-
dence on the chosen observable, but the results
below should hold in the general case of quantum
chaotic or non-integrable systems.
The quantity under study, the time-averaged
fluctuations of the probe observable, is defined
by
δ
2
σ
z
(T ) =
1
T
Z
T
0
dt (hσ
z
(t)i − µ
σ
z
(T ))
2
, (6)
where µ
σ
z
(T ) =
1
T
R
T
0
dthσ
z
(t)i. Assuming
only that the many-body eigenenergies are non-
degenerate, and that also their energy gaps are
non-degenerate, we may express the time fluctu-
ations in terms of matrix elements between eigen-
states of the coupled probe-bath system [25],
δ
2
σ
z
(∞) =
X
µ,ν
µ6=ν
|ρ
µν
|
2
|(σ
z
)
µν
|
2
, (7)
where ρ
µν
= hψ
µ
|ρ|ψ
ν
i, and (σ
z
)
µν
= hψ
µ
|σ
z
|ψ
ν
i.
2.2 Summary of main results of this work
In many practical situations the quantum device
is a non-integrable, quantum chaotic system. We
expect that in such cases a qualitative description
can be obtained by assuming that V is a random
matrix. This assumption leads to a statistical
theory for the many body wave functions in Eqs.
(4, 5) [40],
|ψ
µ
i =
X
α
c
µ
(α)|φ
α
i, (8)
where summations are understood to be taken
from 1 to 2N
B
, the dimension of the Hilbert space
of the device. We will refer to c
µ
(α) as quan-
tum chaotic or random wave functions. We as-
sume that V belongs to the Gaussian Orthogonal
Ensemble (GOE), appropriately scaled by a cou-
pling strength g. In Ref. [40] we showed that this
model can be solved and it allows us to calculate
matrix elements in the interacting basis.
Special care must be taken to translate the
probe initial state defined in Eq. (1) in the ran-
dom matrix formalism. Without loss of general-
ity we can assume that non-interacting eigenfunc-
tions are ordered such that odd (even) values of α
correspond to the probe qubit in state | ↑i (| ↓i).
This leads to the initial condition
ρ(t = 0) =
1
N
B
2N
B
X
α∈odd
|φ
α
ihφ
α
|. (9)
In the following sections we prove that ran-
dom matrix theory leads to a relation between
the time-fluctuations and the typical decay rate
of the probe qubit. This relation is the main re-
sult of this work:
δ
2
σ
z
(∞) = χ(N)Γ
−1
, (10)
where Γ
−1
=
1
∆E
R
∆E
0
dEΓ(E)
−1
is the average
inverse decay rate of the qubit. This average is
an unbiased average over all initially populated
initial state energies, which spans the entire en-
ergy range ∆E of the device due to the initial
infinite temperature bath state. We will see that
this unbiased average originates from the decay
to equilibrium of each state |φ
α
i contributing to
the initial state (9). In Sections 5 and 6 we will
see that this can be related to the decay rate ob-
tained from a fit to the decay to equilibrium of a
local probe observable in our current set up. This
quantity is thus simply an average of the decay
rates experienced by the probe qubit.
The quantity χ(N ), with N the total number
of qubits, depends on the size of the system in
the following way,
χ(N) = C
1
N
B
D(E)
. (11)
D(E) is the average density of states (DOS) of
the system, defined analogously to the average
decay rate above (see also Eq. (61)). N
B
is the
bath Hilbert space dimension, and finally, C is
a constant of order one that does not depend on
the size of the system or coupling strength.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 3
Eq. (10) can be understood as a fluctuation-
dissipation relation, which relates the time-
fluctuations in the steady-state with the decay
rate after a quantum quench. The ratio between
fluctuations and average decay time allows us to
quantify the dimension of the Hilbert space over
which the ergodic quantum dynamics takes place.
In a physical system if V couples the probe qubit
to the whole spectrum of the quantum device, we
expect that the function χ ∝ e
cN
B
, where N
B
is
the number of qubits in the part of the device
that acts as a bath.
Our main result, Eq. (11), can be obtained by
following these steps:
(i) We consider the random matrix model with
an homogeneous coupling matrix V , and calcu-
late δ
2
O
(∞) by extending the formalism developed
in [40, 41] to mixed states (see Section 6.1). To
carry out this calculation we assume a constant
DOS, D(E) = 1/ω
0
, with ω
0
the mean energy
separation.
This calculation allows us to predict an expo-
nential decay for a probe observable of the form:
hO(t)i = (hO(t)i
0
− O
DE
)e
−2Γt
+ O
DE
, (12)
where O
DE
is the diagonal ensemble of the observ-
able O, which we find to be equal to the long-time
average of O as required, and hO(t)i
0
refers to
the free evolution of the observable O under the
Hamiltonian H
0
. This result is derived in Section
5. This is analogous to the result in Reference [41]
for pure-states. We note that the same Equation
has similarly been obtained in Ref. [35], which
allows also for the perturbation matrix V to be
inhomogeneous. In the specific case of interest
here, Eq. (12) becomes,
hσ
z
(t)i = e
−2Γt
, (13)
as hσ
z
(t)i
0
= 1 and hσ
z
i
DE
= 0.
Using this random matrix model we also obtain
a preliminary version of the infinite temperature
fluctuation-dissipation theorem, Eq. (11), with
χ(N) = Cω
0
/N
B
.
(ii) In a realistic system both the DOS and
the coupling strength will depend on the energy.
This leads to energy-dependent qubit decay rates,
Γ(E), and DOS, D(E). Due to the infinite-
temperature initial condition, it is not possible
a priori to approximate those quantities to any
single value, since the quantum dynamics of the
probe qubit result from contributions from all
possible initial states. It is thus necessary to ex-
tend the RMT formalism to allow for variations
in both the DOS and the decay rate in energy
over the width of the initial mixed state.
It is also important to note that when applied
to a realistic system N
B
should be thought of as
the number of eigenstates that may contribute to
the dynamics of a local observable - it is thus a
measure of the effective Hilbert space dimension
in the sense of the number of degrees of freedom
that may be explored from a given initial state.
This measure accounts for effects such as locality
of interactions, and disorder, which will reduce
the Hilbert space available for evolution of the
system. As such, N
B
will in general be bounded
by the total Hilbert space dimension, but is a
more realistic measure of the size of the explored
Hilbert space in the thermalization dynamics of
the device.
In Section 5 we show that the decay rate ob-
served from such an initial state, defined with an
energy dependent coupling strength, is the ther-
mal average decay rate over the energy width de-
fined by the initial state.
We formulate the generalization of the FDT in
Section 6.2, however a brief summary of the ap-
proach is as follows: To account for the energy
variation of the DOS and decay rate we instead
make the assumption that both change slowly in
energy with respect to the energy width of a sin-
gle eigenstate. In this case, one can reformulate
the theory such that the random wave functions
have a smoothly varying width. In effect this is
the statement that a given random wave func-
tion contributes as if it was the eigenstate of a
random matrix model with a constant DOS and
decay rate, as over the width of the wave function
itself these parameters do not change appreciably.
In this case, we see that the average DOS
and decay rate over the initially populated initial
states contributes to the FDT, as in Eq. (10).
This is extended to finite temperature systems,
where instead a thermal average can be seen to
contribute, in Section 6.3.
Finally, we may relate the two decay rate av-
erages, as the thermal average observed from the
decay differs slightly from that appearing in Eq.
(10). These can be seen to be related by a con-
stant factor that does not depend on the coupling
strength, or bath size, and thus scaling informa-
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 4
tion of the Hilbert space dimension may be recov-
ered. This is shown in Section 6.3, and Appendix
E, where the FDT is recast in terms of thermal
averages.
In Section 5 we show that when an exponential
decay of the form Eq. (12) is observed, this in-
dicates that Γ(E) is approximately constant over
the bulk of the initially populated states. In this
case we are able to extract the numerical value
of the Hilbert space dimension directly, as the
averages ove Γ(E) occurring in the FDT, and ob-
served from the decay, are equal.
In summary, we observe the emergence of a
classical fluctuation-dissipation theorem, relating
the time-fluctuations and decay rate of our probe
observable σ
z
. The susceptibility χ(N) in Eq.
(11) can be seen to be related to the Hilbert
space dimension of the bath, N
B
, and thus mea-
surements of the decay rate, Γ, and fluctuations
δ
2
σ
z
(∞), which are both obtainable from the time
evolution, can be exploited to obtain information
on the device Hilbert space dimension, N
B
.
3 Numerical Experiments
Before going into the technical details of our
derivation, we present numerical evidence that
confirms the validity of the random matrix ap-
proach and our main results. In this section, we
show the application to a spin-chain system us-
ing exact diagonalization [42, 43]. From this, we
observe the FDT numerically using a realistic ex-
perimentally observable model.
In Fig. 2, we show the manifestation of Eq.
(10) in a spin-chain system described by the
Hamiltonian H = H
S
+H
B
+H
SB
, where H
S
= 0
is the system Hamiltonian (acting as our probe),
H
B
is our bath Hamiltonian, given by
H
B
=
N
X
j>1
(B
(B)
z
σ
(j)
z
+ B
(B)
x
σ
(j)
x
)+
N−1
X
j>1
[J
z
σ
(j)
z
σ
(j+1)
z
+ J
x
(σ
(j)
+
σ
(j+1)
−
+ σ
(j)
−
σ
(j+1)
+
)],
(14)
which acts on sites with index > 1, which is the
probe index. The probe and bath are coupled by
the interaction Hamiltonian,
H
SB
= J
(SB)
z
σ
(1)
z
σ
(N
m
)
z
+ J
(SB)
x
(σ
(1)
+
σ
(N
m
)
−
+ σ
(1)
−
σ
(N
m
)
+
).
(15)
Figure 2: Observation of the FDT, Eq. (10), for varying
coupling strength, for many bath sizes N − 1 (labelled
on plot). Fits (blue dashed lines) shown for each value
of N, are linear fits to obtain χ(N) for each individual
N value. We see in this case, then, that one does not
require the ability to change the device length in order
to observe our predictions experimentally. Parameters:
B
(B)
Z
= 0.1, B
(B)
x
= 0.3, J
z
= 0.1, J
x
= 1, J
(SB)
z
=
0, β = 0.
Here N
m
is the device site where the probe is
coupled, which we set as 2 throughout. This spin-
chain model may be related to the random matrix
toy model, H = H
0
+ V , via the prescription
H
0
⇔ H
S
+ H
B
, and V ⇔ H
SB
. In particular,
we see that, as χ(N) =
δ
2
O
(∞)
Γ
−1
∝ N
−1
B
, we expect
that if all of the available Hilbert space is being
utilized in the unitary dynamics we will observe
the following scaling:
χ(N) ∝ e
−cN
B
. (16)
This is the relation that we test in Fig. 3.
It is important to note that this exponential
scaling of χ(N), Eq. (16), is expected from not
only the contribution of N
B
, but also from the
average DOS, D(E). This average is often triv-
ially obtained, as for example, for an ensemble of
N two-level systems D(E) =
1
∆E
R
∆E
0
dED(E) =
2
N
∆E
, where ∆E is the range of energies avail-
able E
max
− E
min
(which may itself change
with N), regardless of the microscopic proper-
ties of the DOS. We thus also study the quantity
χ(N)D(E), as this quantity has no dependence
on the DOS, and an observation of the exponen-
tial scaling in system size is confirmation that,
indeed, N
B
∝ e
cN
. This is shown in Fig. 4,
where we observe an exponential scaling of the
Hilbert space dimension, with c ≈ 0.62, compared
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 5
Figure 3: Observation of the FDT, Eq. (10), for vary-
ing bath size N − 1, for temperatures β. Fit (yellow
dashed line) is performed to the function ae
N
N
0
for the
infinite temperature case, β = 0, and thus confirms the
exponential scaling of χ(N). In this case, we observe an
exponential scaling for all temperatures, as the average
DOS also scales exponentially with N. Note that here
we have Γ ∼ 0.2 so the high temperature limit is defined
by approximately β 5. Parameters: as Fig. 2 with
J
(SB)
x
= 0.5.
to ln(2) ∼ 0.69 if the entire Hilbert space were
explored in the dynamics.
We further observe in Figs. 3 and 4, that
the FDT similarly applies at finite temperatures
β =
1
k
B
T
> 0. The extension of our theoretical
approach to this case is discussed below, with ad-
ditional details given in Appendix E. Indeed, we
can show that for high temperatures, such that
β
−1
Γ, we obtain an FDT of the same form
as Eq. (10), by employing a high energy cut-off
(ρ
B
)
αα
∼ e
−βE
α
to the bath state occupation.
For finite temperatures we show below that the
FDT depends on the partition function Z
β
itself,
rather than the Hilbert space dimension. Indeed,
one can see that in the infinite temperature limit
Z
0
= lim
β→0
P
α
e
−βE
α
δ
α,odd
= N
B
.
Finally, we see that when the observable is
found to decay exponentially to its equilibrium
value, this indicates that the decay rate is ap-
proximately constant over the bulk of the ini-
tially occupied states. This is shown in Section
5. Exploiting this observation, we are able to di-
rectly obtain the Hilbert space dimension, as for
an infinite temperature initial bath state the av-
erage Γ
−1
is equal to the measured decay rate.
The bath Hilbert space dimension N
B
, as calcu-
lated from Eq. (10), is plotted for varying de-
Figure 4: As in Fig. 3, however accounting for the expo-
nential scaling of the average DOS with device size. Fit
(yellow dashed line) is performed to the function ae
N
N
0
for the infinite temperature case, β = 0, and thus con-
firms the exponential scaling of χ(N)D(E) ∼
1
N
B
. In
this case, we observe an exponential scaling for only
high temperatures, and confirm that for low tempera-
tures χ(N)D(E) is independent of N. Note that here
we have Γ ∼ 0.2 so the high temperature limit is defined
by approximately β 5. Parameters: as Fig. 2 with
J
(SB)
x
= 0.5.
vice sizes in Fig. 5. Here we observe that N
B
indeed increases exponentially with systems size,
yet is somewhat smaller than its maximum possi-
ble value 2
N−1
, which is expected to the locality
of interactions within the chain.
We again note that in Fig. 5 the measurement
of N
B
is a measurement of the explored Hilbert
space, or the total number of eigenstates that con-
tribute to the evolution of the initial state. Thus,
for a maximally connected device this would be
2
N−1
, whereas locality of interactions in this case
restricts some areas of the Hilbert space.
We note that for models where Γ(E) is not ap-
proximately constant in energy, which would be
marked by a deviation from the exponential de-
cay in Eq. (12), one would still have access to
Figs. 2-4, and thus scaling information of the
Hilbert space dimension, yet the numerical value
N
B
would be obscured. This is explained in more
detail in Section 5, where we see how information
on Γ(E) is extracted from the decay, and in Sec-
tion 6 and Appendix E, where we see how this
relates to the observed FDT. The key detail is
that when the FDT is expressed in terms of the
same value measured from the decay it differs by
a constant, thereby leaving the scaling of N
B
with
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 6
Figure 5: Hilbert space dimension of the bath N
B
for
varying device size N inferred from the FDT, Eq. (10)
(blue dots), C =
3
8π
, which can be easily obtained from
the observable (see Section 6). For a fully connected
device the maximum possible N
B
is N
B
= 2
N−1
, the
region below this limit is shaded. We observe that the
measured Hilbert space dimension indeed scales expo-
nentially with system size (dashed line is exponential fit),
though is somewhat smaller than the maximum dimen-
sion of a fully connected device. Parameters as in Fig.
2, with J
(SB)
x
= 0.5 and β = 0.
system size or decay rate accessible, yet obscuring
the numerical value.
We note that in Ref. [41], the current authors
obtained a FDT for pure states, which can be
seen to be recovered in the low temperature limit,
β
−1
Γ, for which χ(N) does not depend explic-
itly on the Hilbert space dimension N
B
. This can
also be analytically seen to be the same as the low
temperature limit of our treatment below, which
indicates that there is a smooth transition be-
tween these two cases. This is indeed observed in
the numerics of Figs 3 and 4.
4 Model
4.1 RMT Approach
Our approach relies on the calculation of correla-
tion functions from a statistical theory of random
wave functions c
µ
(α). Here we summarize the es-
sential ingredients to our model, and give details
on the calculations in the sections below.
Our theory, developed in Ref. [40] by ex-
tending Deutch’s RMT model [44–46], can be
used to obtain arbitrary correlation functions
hc
µ
(α)c
ν
(α) · · · i
V
, where h· · · i
V
denotes the en-
semble average over an ensemble of random ma-
trix perturbations, V , for a N × N Hamiltonian
of the form (3), with (H
0
)
αβ
= αω
0
δ
αβ
, with
ω
0
=
1
N
and V a random matrix selected from the
GOE, with hV
2
αβ
i
V
=
(1+δ
αβ
)g
2
N
. In practice, those
correlations allow us to calculate any dynamical
quantity of interest within the RMT formalism.
To illustrate the use of such correlation func-
tions we briefly consider the simple example of
the diagonal observable matrix elements O
µµ
.
These can be written as
O
µµ
=
X
αβ
c
µ
(α)c
µ
(β)O
αβ
. (17)
Now, using the self-averaging property of random
matrices, which we prove for this model in Ap-
pendix C, we see that
O
µµ
=
X
αβ
hc
µ
(α)c
µ
(β)i
V
O
αβ
. (18)
Note that only the random wave functions c
µ
(α)
remain inside the ensemble average, as all other
factors are in the non-interacting basis and thus
do not depend on V . We thus observe that the di-
agonal observable matrix elements depend on the
correlation function hc
µ
(α)c
µ
(β)i
V
(note that we
will see how to deal with the summation and non-
interacting observable elements in Section 4.3 be-
low).
More generally, when calculating more compli-
cated quantities, we have that there is also a non-
trivial contribution of a four-point correlation
function of the form hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
)i
which is given by (for µ 6= ν),
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 7
hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
)i
V
= Λ(µ, α)Λ(ν, β)δ
αα
0
δ
ββ
0
−
Λ(µ, α)Λ(ν, β)Λ(µ, α
0
)Λ(ν, β
0
)
Λ
(2)
(µ, ν)
(δ
αβ
δ
α
0
β
0
+ δ
αβ
0
δ
βα
0
),
(19)
where Λ(µ, α) is defined as
hc
µ
(α)c
ν
(β)i
V
=Λ(µ, α)δ
αβ
δ
µν
(20a)
Λ(µ, α) :=
ω
0
Γ/π
(E
α
−E
µ
)
2
+Γ
2
(20b)
with Γ =
πg
2
N ω
0
. Λ
(2)
(µ, ν) is defined similarly to
Eq. (20b), with Γ → 2Γ. The Lorentzian form
of Eq. (20) is found for a homogeneous perturba-
tion V [40, 45], selected from the GOE. We ob-
tain that the four-point correlation function, Eq.
(19), can be described in terms of product of two-
point correlators Λ(µ, α), as if the random wave
function distribution was purely Gaussian, plus a
correction term originating from the effective in-
teraction due to the mutual orthogonality of ran-
dom wave functions. We will see below that an
approach in terms of Gaussian and non-Gaussian
contractions can be formulated to describe more
general correlation functions.
We note that our theory can be extended to ac-
count for different forms of the quantum chaotic
wave function Λ(µ, α). This may appear, for ex-
ample, for non-homogeneous V . In this case, the
form of the function Λ would change, however the
algebraic structure of our theory would remain.
In order to evaluate Eq. (7), we thus use
Eq. (8) to write δ
2
O
(∞) in terms of the ran-
dom wave functions c
µ
(α), and non-interacting
matrix elements ρ
αβ
and (σ
2
z
)
αβ
. We then use
the self-averaging property of random matrices
(which we discuss in A, and prove in Appendix
C), and obtain the relevant correlation functions
hc
µ
(α)c
ν
(β) · · · i
V
.
4.2 Computing Correlation Functions
As we have seen, it is important to have a sys-
tematic approach to obtaining correlation func-
tions for this model. This is a non-trivial task
as the random wave functions of our theory are
not Gaussian independent variables, but include
an effective interaction due to the orthogonality
condition hψ
µ
|ψ
ν
i = δ
µν
[40]. Our approach to
the statistical theory of random wave functions is
summarized in Appendix A.
Below, we present such a systematic approach
to obtaining arbitrary correlation functions in
terms of contractions representing the Gaussian
and non-Gaussian terms in the four-point corre-
lator (which is the largest non-factorizable corre-
lation function of our theory).
The four-point correlation function of Eq.
(19) may be understood in terms of the con-
tractions of non-interacting indices, indeed
it can be seen to be the sum of a Gaussian
contraction hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
)i
V
⇒
hc
2
µ
(α)i
V
hc
2
ν
(β)i
V
δ
αα
0
δ
ββ
0
=
Λ(µ, α)Λ(ν, β)δ
αα
0
δ
ββ
0
and non-Gaussian con-
tractions, given by
hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
)i
V
⇒ L
αβα
0
β
0
µν
δ
αβ
δ
α
0
β
0
hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
)i
V
⇒ L
αβα
0
β
0
µν
δ
αβ
0
δ
α
0
β
,
(21)
where
L
αβα
0
β
0
µν
:=
Λ(µ, α)Λ(ν, β)Λ(µ, α
0
)Λ(ν, β
0
)
Λ
(2)
(µ, ν)
. (22)
We reserve the double line contraction notation
of Eq. (21) for the non-Gaussian case. Note that
these must occur in pairs of contractions between
different interacting indices µ 6= ν.
Now, we can see from Eq. (21) that each con-
tracted pair of indices contributes a Kronecker-
δ symbol, and thus, when the correlation func-
tion is summed over its non-interacting indices,
the number of summations is reduced. We see
that as each Λ contributes a factor on the or-
der O(
ω
0
Γ
), and a summation on the order O(
Γ
ω
0
),
a reduced summation will act to render a term
negligible in comparison to a term with no such
reduction. Further, we see that the contribu-
tion of the non-Gaussian term Eq. (22) is of or-
der O(
ω
3
0
Γ
3
), whereas that of the Gaussian term is
∼ Λ
2
, and thus O(
ω
2
0
Γ
2
), and as such, one can see
that for the non-Gaussian contractions to con-
tribute, they must be acted on my an extra sum-
mation. Indeed, one can see that this occurs for
one of the two non-Gaussian terms when one has
repeated summations, i.e. α
0
, β
0
→ α, β in Eq.
(21).
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 8
For further details we refer the reader to Ref.
[41], and the calculations in Appendices B and
C. Here we have seen the key intuition, how-
ever: that repeated indices in correlation func-
tions leads to the dominant contribution of con-
tractions that would otherwise have contracted
the pair of equal indices.
4.3 Assumptions on Observables
After obtaining the relevant correlation function,
one needs to perform the summations over re-
maining indices. See, for example, the simple
case of Eq. (18) above. To perform the sum-
mations, certain assumptions on the form of ob-
servable matrix elements in the non-interacting
basis must be made. We note that in this basis,
local system observables are usually of a known
form.
Our theory relies on assumptions that we ex-
pect to be satisfied for such local observables.
The key assumption is related to the behaviour of
matrix elements, which must have a well defined
average, that does not vary pathologically with
energy. A more formal definition of our assump-
tion can be written in terms of the function Λ,
as,
X
α
Λ(µ, α)Λ(ν, α)O
αα
= [O
αα
]
µ
Λ
(2)
(µ, ν), (23)
with E
µ
:=
E
µ
+E
ν
2
, and
[O
αα
]
µ
=
X
α
Λ(µ, α)O
αα
. (24)
We will see that this assumption will be necessary
in order to compute summations over the non-
interacting indices. In this section we explain in
more detail the requirements on the form of O
αα
for Eq. (23) to be valid, as well as the physical
interpretation of the assumption.
The essential assumption here, which we label
smoothness of O
αα
, as in Ref. [41], is that the mi-
crocanonical average [O
αα
]
µ
changes slowly over
the width Γ of the function Λ(µ, α)Λ(ν, α). We
showed in Ref. [41] that this is the case under the
assumptions,
Γ
ω
0
1,
Γ
2
d
2
dE
2
µ
[O
α,α
]
µ
1, (25)
which thus leads us to two reasonable conditions,
1. There are many states in the energy width Γ
2. The microcanonical average changes slowly
over the width Γ.
We note that the latter condition, combined with
the fact that the microcanonical average and time
average are equal (which is shown below), is
equivalent to the statement that the time-average
of the observable is not sensitive to the particular
initial state (microstate), rather, it’s macroscopic
energy. In fact, one can see that the conditions
(25) are precisely those required in order to de-
fine a microcanonical average that does not vary
pathologically with small changes in the energy
window. In this sense, this assumption is the min-
imal assumption one would expect to require for
thermalization to a microcanonical average that
changes smoothly with initial state energy to oc-
cur.
We further note that in the consideration of
time evolution below, we will consider more gen-
eral observables that are not necessarily diagonal
in the non-interacting basis, but fulfil a sparsity
condition. This can be written as
P
αβ
O
αβ
=
P
α
P
N
O
n
O
α,α+n
δ
β,α+n
, where for a given ob-
servable there is a non-extensive number N
O
of
groups of non-zero matrix elements at given en-
ergy widths, such that after the course grain-
ing procedure the observable matrix elements are
non-zero for energy gaps E
α
−E
β
that are the pos-
sible energy gaps of H
S
. This form can be seen
[41] to be reasonable for local observables. We
note that it is of course possible to find observ-
ables that do not fulfil this assumption, although
it is easily seen to be true for e.g. local Pauli
operator observables. We will see below that our
treatment of time evolution could potentially also
capture a wider range of observables as well, if
the form in the non-interacting basis is known.
In the following we refer to observables fulfilling
the above assumptions as ‘generic’ observables.
5 Equilibration Dynamics
In this section we present a description of the
time dependence of ‘generic’ observables as de-
fined above, from an arbitrary initial condition
ρ(0) =
X
αβ
w
αβ
|φ
α
ihφ
β
|. (26)
This calculation may be performed by exploiting
the methods outlined in Section 4. The general
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 9
approach may be summarized in three steps: i)
Writing the observable time dependence in terms
of parameters in the non-interacting basis, ii)
computing the relevant correlation functions (see
Section 4.2), and iii) performing summations us-
ing the assumptions on observables (see Section
4.3).
Proceeding as such, we write the time depen-
dent density operator in the form,
ρ(t) =
X
αβ
X
µν
w
αβ
c
µ
(α)c
ν
(β)e
−i(E
µ
−E
ν
)t
|ψ
µ
ihψ
ν
|,
(27)
which may be used to obtain the time evolved
observable expectation value by hO(t)i =
Tr(ρ(t)O). By taking the trace over the inter-
acting basis {|ψ
µ
i}, we thus obtain
hO(t)i =
X
µ
0
hψ
µ
0
|
X
αβµν
w
αβ
c
µ
(α)c
ν
(β)
× e
−i(E
µ
−E
ν
)t
|ψ
µ
ihψ
ν
|O|ψ
µ
0
i.
(28)
Noting the so-called diagonal ensemble contribu-
tion is defined by,
O
DE
=
X
αµ
w
αα
c
2
µ
(α)O
µµ
, (29)
which can be seen to be equal to the long-time
average value of the observable
hhO(t)ii
τ
:= lim
τ→∞
1
τ
Z
τ
0
dthO(t)i = O
DE
, (30)
assuming no degenerate energy levels, we thus de-
fine
∆O(t) : = hO(t)i − O
DE
=
X
αβ
X
µν
µ6=ν
w
αβ
e
−i(E
µ
−E
ν
)t
O
µν
× c
µ
(α)c
ν
(β).
(31)
Using that O
µν
=
P
αβ
c
µ
(α)c
ν
(β)O
αβ
, we have
∆O(t) =
X
αβα
0
β
0
X
µν
µ6=ν
w
αβ
e
−i(E
µ
−E
ν
)t
O
α
0
β
0
× c
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
).
(32)
We see that, using the self-averaging property,
this depends on the four-point correlation func-
tion given by Eq. (19), such that
∆O(t) =
X
µν
µ6=ν
e
−i(E
µ
−E
ν
)t
×
"
X
αβ
w
αβ
O
αβ
Λ(µ, α)Λ(ν, β)
−
X
αβ
w
αα
O
ββ
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
−
X
αβ
w
αβ
O
αβ
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
#
.
(33)
The third term can be shown to be negligible,
proof of which is given in Appendix D. We note
that it is this bound that requires the sparsity
assumption above. We now use the smoothness
assumption, exploiting Eqs. (23) and (24), we
obtain,
∆O(t) = hO(t)i
0
e
−2Γt
−
X
µν
µ6=ν
X
α
[O
αα
]
µ
w
αα
e
−i(E
µ
−E
ν
)t
Λ(µ, α)Λ(ν, α),
(34)
where to obtain the first term one may note that
Λ(µ, α) = Λ(µ − α) = Λ(α − µ), and make the
change of variables µ(ν) → µ(ν) − α(β) to per-
form the integrals over the new variables. Here
hO(t)i
0
:=
P
αβ
w
αβ
O
αβ
e
−i(E
α
−E
β
)t
is the free
evolution of the observable under the Hamilto-
nian H
0
. The second term may be re-expressed
by defining ˜µ = µ − α, to obtain,
∆O(t) = hO(t)i
0
e
−2Γt
−
X
˜µ˜ν
˜µ6=˜ν
X
α
[O
αα
]
m
w
αα
e
−i(E
˜µ
−E
˜µ
)t
Λ(˜µ)Λ(˜ν),
(35)
where E
m
:=
E
˜µ
+E
˜ν
2
+ E
α
. Noting, then, that as
Λ(˜µ) is peaked around E
˜µ
= 0, and that [O
αα
]
α
changes slowly over the width Γ of the function Λ,
we can make the replacement [O
αα
]
m
→ [O
αα
]
α
.
This allows the summations over ˜µ, ˜ν to be per-
formed, which become Fourier transforms of the
Lorentzian functions Λ in the continuum limit
P
µ
→
R
dE
µ
ω
0
. We thus find,
∆O(t) = hO(t)i
0
e
−2Γt
−
X
α
[O
αα
]
α
w
αα
e
−2Γt
.
(36)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 10
Noting that at t = 0 we by definition
have ∆O(0) := hO(0)i − O
DE
= hO(0)i
0
−
P
α
[O
αα
]
α
w
αα
, we obtain that
P
α
[O
αα
]
α
w
αα
=
O
DE
. Noting Eq. (30), we see that the equality
of the time and microcanonical averages is de-
rived from our RMT approach. Thus, using the
definition in Eq. (31), we obtain
O(t) = (hO(t)i
0
− O
DE
)e
−2Γt
+ O
DE
. (37)
This is the same as that obtained in Reference
[41] for pure-states.
The approach outlined above is valid assuming
that the decay rate Γ is constant in energy. In
fact, for the system under consideration, we have
to allow Γ to change with the initial state energy,
such that Γ → Γ
α
(note that here the change in
DOS does not affect the calculation). Accounting
for this, rather than Eq. (34), we obtain
∆O(t) =
X
αβ
w
αβ
O
αβ
e
−i(E
α
−E
β
)t
e
−(Γ
α
+Γ
β
)t
−
X
α
[O
αα
]
α
w
αα
e
−2Γ
α
t
.
(38)
Now, for our system we have that, as H
S
= 0,
the microcanonical average [O
αα
]
µ
= 0 for all µ
in the bulk of the spectrum. Also, using that
for our proposed experimental protocol, both the
initial state and observable are diagonal in the
non-interacting basis, we have
∆O(t) =
X
α
w
α
O
αα
e
−2Γ
α
t
. (39)
We have that, as the initial state hO(0)i =
P
α
w
α
O
αα
= 1 for all initial device states, O
αα
=
1 for all non-zero w
α
, and thus
∆O(t) =
X
α
w
α
e
−2Γ
α
t
. (40)
We thus wish to obtain the value Γ that will be
obtained when measuring the decay of an observ-
able. To find this, one may simply consider the
time integration of the evolution obtained above
from the initial state,
ρ
α
=
1
Z
β
e
−βE
α
δ
α,odd
, (41)
describing our probe-bath model, with an initial
finite temperature bath state at inverse tempera-
ture β. The time integration is then,
lim
τ→∞
Z
τ
0
dt∆O(t)
= lim
τ→∞
Z
τ
0
dt
X
α
w
α
e
−2Γ
α
t
= −
1
2
X
α
Γ
−1
α
w
α
= −
1
2Z
β
X
α
Γ
−1
α
e
−βE
α
δ
α,odd
:= −
1
2
hhΓ(E)
−1
ii
β
,
(42)
where we have used in the second line that O
DE
=
0, and defined the thermal average hh· · · ii
β
at
inverse temperature β, and we have defined Γ
α
as
the decay rate of the initial state |φ
α
i. We thus
see that it is the thermal average of the inverse
decay rate that is measured by a fit to the time
dependence of an observable.
The integral form of the thermal average of a
function A(E), is given by,
hhA(E)ii
β
:=
1
Z
0
β
Z
∆E
0
dED(E)e
−β(E−E
0
)
A(E),
(43)
with Z
0
β
:=
R
∆E
0
dED(E)e
−β(E−E
0
)
. We will see
below that the FDT will be initially expressed in
terms of a different average over the values Γ(E).
This difference is resolved in Appendix E, where
we re-express our FDT in terms of the thermal
average above. We show that the form differs
only by a constant that is independent of N and
the coupling strength, and thus the scaling with
Hilbert space dimension remains the same, and
this difference is not important for our applica-
tion.
Finally, we note that when Γ
α
≈ Γ is ap-
proximately constant across the bulk of the ini-
tially populated initial states, the thermal aver-
age above is approximately equal to the unbiased
average appearing in (10). We also see that in
this case, from Eq. (40), one expects to observe
an exponential decay at the rate Γ, as in Eq. (13).
Indeed, this is what we observe in our numerical
example in Section 3 above, and thus we are able
to recover the Hilbert space dimension directly
in Fig. 5. If a non-exponential decay is observed,
then the average hhΓ(E)
−1
ii
β
is obtainable via in-
tegration as above, and the scaling of the Hilbert
space dimension is still obtainable as in Fig. 4.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 11
6 Fluctuation-Dissipation Theorem
6.1 Derivation from RMT
Here we perform the full derivation of the FDT
for the random matrix model described above.
We initially focus on the case of a diagonal ini-
tial bath state ρ
αβ
= w
α
δ
αβ
. We then restrict
the treatment to the specific protocol outlined in
Section 2, where the initial state is the product of
a single probe qubit in a pure state, and a bath
in an infinite temperature state, see Eq. (9). We
will follow a very similar steps as those outlined
in the previous section, however we will see here
that the correlation function calculation is some-
what more complicated.
The RMT model here is limited to the case of
constant decay rate and DOS, we will thus extend
the treatment to more realistic cases in the next
section.
We are interested in the calculation of the long-
time fluctuations, defined by the diagonal ensem-
ble result,
δ
2
σ
z
(∞) =
X
µ,ν
µ6=ν
|ρ
µν
|
2
|(σ
z
)
µν
|
2
.
(44)
We begin be writing the initial density operator
matrix elements as,
ρ
µν
=
X
α
w
α
c
µ
(α)c
ν
(α), (45)
then using Eqs. (44), and that |ψ
µ
i =
P
α
c
µ
(α)|φ
α
i, we may write the time fluctuations
as,
δ
2
O
(∞) =
X
µν
µ6=ν
X
αβ
w
α
w
β
c
µ
(α)c
ν
(α)c
µ
(β)c
ν
(β)
×
X
α
0
β
0
O
α
0
α
0
O
β
0
β
0
c
µ
(α
0
)c
ν
(α
0
)c
µ
(β
0
)c
ν
(β
0
),
(46)
where coefficients of the initial state are labelled
as unprimed indices α, β, and coefficients of the
observable are labelled by primed indices.
Using the self-averaging property of random
matrices (see Appendices A and C), we may re-
place the product of coefficients c
µ
(α)c
ν
(α) · · ·
by their ensemble average hc
µ
(α)c
ν
(α) · · · i
V
; the
above expression may then be written in terms
of a sum over 8-point correlation functions,
weighted by the initial state and observable co-
efficients w
α
and O
αα
:
δ
2
O
(∞) =
X
µν
µ6=ν
X
αβα
0
β
0
w
α
w
β
O
α
0
α
0
O
β
0
β
0
hc
µ
(α)c
ν
(α)c
µ
(β)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
)c
µ
(β
0
)c
ν
(β
0
)i
V
. (47)
Now, using the method of contractions outlined
in Section 4.2, we see that this 8-point correlation
function may be split up into to four-point corre-
lation functions, each consisting of both Gaussian
and non-Gaussian contractions. These are com-
puted explicitly in Appendix B, in which we see
that there are three dominating contributions to
the fluctuations, given by,
δ
2
G
(∞) =
X
µν
µ6=ν
[w
2
α
]
µ
[O
2
αα
]
µ
+ 2[w
α
O
αα
]
2
µ
Λ
(2)
(µ, ν)
2
,
(48)
δ
2
NG
(∞) = 3
X
µν
µ6=ν
[w
α
]
2
µ
[O
αα
]
2
µ
Λ
(2)
(µ, ν)
2
, (49)
and,
δ
2
M
(∞) = −
X
µν
µ6=ν
Λ
(2)
(µ, ν)
2
[w
2
α
]
µ
[O
αα
]
2
µ
+ 4[w
α
]
µ
[O
αα
]
µ
[w
α
O
αα
]
µ
+ [w
α
]
2
µ
[O
2
αα
]
µ
.
(50)
These three terms can be seen as the contri-
butions to the 8-point correlation function aris-
ing due to products of Gaussian, non-Gaussian,
and mixed Gaussian and non-Gaussian 4-point
correlation functions respectively. In the above,
in order to perform the summations over non-
interacting indices in Eq. (47) we define course
grained averages of observable elements O
αα
as
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 12
in Eq. (23), as well as the mixed averages,
X
α
Λ(µ, α)Λ(ν, α)w
α
= [w
α
]
µ
Λ
(2)
(µ, ν)
X
α
Λ(µ, α)Λ(ν, α)w
α
O
αα
= [w
α
O
αα
]
µ
Λ
(2)
(µ, ν).
(51)
We thus define W
µ
by
δ
2
O
(∞) = δ
2
G
(∞) + δ
2
NG
(∞) + δ
2
M
(∞)
=
X
µν
µ6=ν
W
µ
Λ
(2)
(µ, ν)
2
,
(52)
with,
W
µ
= [w
2
α
]
µ
[O
2
αα
]
µ
+ 2[w
α
O
αα
]
2
µ
+ 3[w
α
]
2
µ
[O
αα
]
2
µ
− [w
2
α
]
µ
[O
αα
]
2
µ
− 4[w
α
]
µ
[O
αα
]
µ
[w
α
O
αα
]
µ
− [w
α
]
2
µ
[O
2
αα
]
µ
.
(53)
We now take our bath to be in an initial infi-
nite temperature state, such that [w
α
]
µ
= [w
α
] =
1
2N
B
, and [w
2
α
]
µ
= [w
2
α
] =
1
2N
2
B
. As such, W
µ
= W
is in fact energy independent, as the probe Hamil-
tonian H
S
= 0, so microcanonical averages of
probe observables are also energy independent.
Now, as [w
2
α
] = 2[w
α
]
2
, all terms in W are ∝ [w
2
α
],
we define,
W
O
=
W
[w
2
α
]
= [O
2
αα
] + O
2
↑
+
3
2
[O
αα
]
2
− [O
αα
]
2
− 2[O
αα
]O
↑
−
1
2
[O
2
αα
],
(54)
where O
↑
= h↑ |O| ↑i, and we have used that
[w
α
O
αα
] = O
↑
[w
α
]. We see that W
O
is a con-
stant of the order of unity that depends only on
the observable (e.g. W
O
= 3/2 for O = σ
z
). Fi-
nally, taking the thermodynamic limit, such that
P
µν
→
R
∆E
0
R
∆E
0
dE
µ
dE
ν
ω
2
0
(not that the diagonal
terms in the summation can be seen to be negli-
gible, as they contribute to a higher order in
ω
0
Γ
),
we have
δ
2
O
(∞) =
W
O
2N
2
B
Z
∆E
0
Z
∆E
0
dE
µ
dE
ν
ω
2
0
Λ
(2)
(µ, ν)
2
,
(55)
which may be evaluated using,
Z
∆E
0
Z
∆E
0
dE
µ
dE
ν
ω
2
0
Λ
(2)
(µ, ν)
2
=
Z
∆E
0
dE
µ
arccot
2Γ
∆E
2πΓ
≈
∆E
4πΓ
,
(56)
where in the last line we have used that ∆E Γ,
such that arccot
2Γ
∆E
≈
π
2
. We then obtain,
δ
2
O
(∞) =
W
O
ω
0
4πN
B
Γ
,
(57)
where we have used that ∆E := 2N
B
ω
0
.
We can see, then, that Eq. (57) is of the form of
our main result, Eq. (3) of the main text, where
C =
W
O
4π
. What follows is to generalize this rela-
tion, allowing the DOS and Γ to vary in energy,
and for finite temperatures.
6.2 Extension to Realistic Systems
The key issue with directly applying the RMT
results to realistic models is that in general the
DOS, and decay rate, are energy dependent, and
thus change over the width of the initial state
distribution (this is especially important for the
high/infinite temperatures considered here). In
order to account for this, we must then go back
to the evaluation of the integrals over energy,
in Eq. (55), and substitute Γ → Γ(E), and
ω
−1
0
→ D(E). This is justified under the assump-
tion that neither Γ(E), nor D(E), vary apprecia-
bly over the width Γ. i.e.
Γ(E)−Γ(E+Γ(E))
Γ(E)
1,
and
D(E)−D(E+Γ(E))
D(E)
1.
We see then, that the integral in Eq. (55) is
now
δ
2
O
(∞) =
W
O
2N
2
B
Z
∆E
0
Z
∆E
0
dE
µ
dE
ν
×
D(E
µ
)D(E
ν
)
D(E)
2
2Γ(E)/π
(E
µ
− E
ν
)
2
+ 4Γ(E)
2
!
2
,
(58)
where we have used that,
Λ
(2)
(µ, ν) =
1
D(E)
2Γ(E)/π
(E
µ
− E
ν
)
2
+ 4Γ(E)
2
, (59)
with E =
E
µ
+E
ν
2
. This can be seen to be valid as
long as the above conditions on D(E) and Γ(E)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 13
hold, that is, as long as they change sufficiently
slowly over the energy width Γ. The contribution
of each eigenstate Λ(µ, α) to the fluctuations at
a given energy is then that of a local (in energy)
random matrix model, with a constant DOS and
decay rate. D(E) and Γ(E) can then be allowed
to change over an energy much wider than a the
width of Λ(µ, α), as over such energy widths the
contributions of relevant eigenstates are indepen-
dent.
Now, we further define ω = E
µ
−E
ν
, and make
the change of variables E
µ
, E
ν
→ ω, E, and thus
obtain
δ
2
O
(∞) =
W
O
2N
2
B
Z
∆E
0
Z
∆E
−∆E
dEdω
×
D(E +
ω
2
)D(E −
ω
2
)
D(E)
2
2Γ(E)/π
ω
2
+ 4Γ(E)
2
2
≈
W
O
2N
2
B
Z
∆E
0
dE
1
4πΓ(E)
,
(60)
where in the second line we have assumed that
D(E) and Γ(E) is approximately constant over
the width Γ. Now, we define the unbiased average
of a function A(E) as,
A(E) =
1
∆E
Z
∆E
0
dEA(E), (61)
(not to be confused with the average [· · · ]
µ
above)
and see that, noting ∆E =
2N
B
D(E)
,
δ
2
O
(∞) =
W
O
4πN
B
D(E)
Γ(E)
−1
= C
1
N
B
D(E)
Γ(E)
−1
,
(62)
where C =
W
O
4π
depends only on the choice of
observable. We note that for the random ma-
trix model, as the DOS and Γ(E) are both con-
stant in energy, the average Γ(E)
−1
is equal to the
thermal average hhΓ(E)
−1
ii
β=0
obtained from a
fit to the decay of an observable (see Section 5
above). In the case above, however, where the
DOS and Γ(E) change in energy, the unbiased
average decay rate is not necessarily the same as
that obtained from a fit to the decay. We may
fix this problem directly, as we do in the last sec-
tion, where we see that the unbiased thermal av-
erages may be replaced by regular thermal aver-
ages weighted by the DOS at the expense of a
constant that depends on the functional form of
D(E) and Γ(E) (but importantly, not on N, or
the coupling strength). We can also see, that if
Γ(E) is approximately constant over the width of
the DOS, which is often the case in such systems
(in fact from Eq. (40) this can be seen to be the
case if an exponential decay of the observable is
observed), then the biased and unbiased thermal
averages of Γ(E)
−1
are approximately equal for
β → 0, and Eq. (62) may be directly experimen-
tally confirmed as in Fig. 5.
6.3 Finite Temperature FDT
In this section we extend the above approach to
finite temperature initial bath states, where the
initial state is described by
ρ(t = 0) =
2N
B
X
α
w
α
|φ
α
ihφ
α
|, (63)
where the joint probe-bath Hamiltonian eigenba-
sis is built by ordering product states such that
|φ
α
i = | ↑i|φ
B,
α+1
2
i (α odd) , |φ
α
i = | ↓i|φ
B,
α
2
i
(α even). In this case, we have
w
α
=
1
Z
β
e
−βE
α
if α ∈ odd
0 otherwise
, (64)
where Z
β
=
P
α
e
−βE
α
δ
α,odd
, when the bath is
initially a finite temperature state at inverse tem-
perature β =
1
k
B
T
, and the probe qubit is initially
in state | ↑i. We thus obtain for the microcanon-
ical averages of w
α
, assuming that β
−1
Γ,
[w
α
]
µ
=
1
2Z
β
e
−βE
µ
(65)
and
[w
2
α
]
µ
=
1
2Z
2
β
e
−2βE
µ
, (66)
such that [w
2
α
]
µ
= 2[w
2
α
]
µ
. Now, our most gen-
eral form for the long-time fluctuations (which
assumes only the ability to define the required
microcanonical averages that vary smoothly over
a width Γ) is
δ
2
O
(∞) =
X
µν
µ6=ν
W
µ
Λ
(2)
(µ, ν)
2
,
(67)
where E
µ
:=
E
µ
+E
ν
2
, and W
µ
is written in Eq.
(53). Indeed, noting that the mixed average
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 14
[w
α
O
αα
]
µ
= 2[w
α
]
µ
O
↑
, where O
↑
:= h↑ |O| ↑i,
and that as each term in W
µ
is ∝ [w
α
]
2
µ
(see Eq.
(53)), we may define
W
µ
= W
O
[w
2
α
]
µ
=
W
O
2Z
2
β
e
−2βE
µ
, (68)
so
δ
2
O
(∞) =
W
O
2Z
2
β
X
µν
µ6=ν
e
−2βE
µ
Λ
(2)
(µ, ν)
2
.
(69)
This may be evaluated, including a variable DOS
D(E), as in the main text for the infinite temper-
ature case, via
δ
2
O
(∞) =
W
0
2Z
2
β
Z
∆E
0
Z
∆E
−∆E
dEdωe
−2Eβ
×
D(E +
ω
2
)D(E −
ω
2
)
D(E)
2
2Γ(E)/π
ω
2
+ 4Γ(E)
2
2
≈
W
O
4Z
2
β
Z
∆E
0
dEe
−2βE
1
4πΓ(E)
,
(70)
where in the second line we have made the change
of variables E
µ
, E
ν
→ E, ω with E =
E
µ
+E
ν
2
and
ω = E
µ
− E
ν
, and used that ∆E Γ as in the
main text. We now define the unbiased thermal
average of the function A(E) as,
hA(E)i
β
:=
1
∆E
0
(β)
Z
∆E
0
dEA(E)e
−βE
, (71)
where ∆E
0
(β) =
R
∆E
0
dEe
−βE
. Now, we have
δ
2
O
(∞) =
W
O
∆E
0
(β)
8πZ
2
β
hΓ(E)
−1
i
2β
. (72)
Noting, then, that lim
β→0
Z
β
= N
B
, and
lim
β→0
∆E
0
(E) = ∆E =
2N
B
D(E)
, we recover the
infinite temperature case as required,
δ
2
O
(∞) =
W
O
4πN
B
D(E)
Γ(E)
−1
. (73)
We note that, unlike in the RMT case above, the
average hΓ(E)
−1
i
β
is not equal to the thermal av-
erage hhΓ(E)
−1
ii
β
, which is that obtained from a
fit to the decay. In Appendix E we show how the
FDT may be defined in terms of this thermal av-
erage. Importantly for our proposed application,
we obtain that the FDT in this form is related
simply by a constant C
0
β
, defined in Appendix E,
that does not depend on the the size of the device
or coupling strength (within the weak coupling
regime). For infinite temperatures we thus have
δ
2
O
(∞) = C
0
0
W
O
4πN
B
hhD(E)ii
0
hhΓ(E)
−1
ii
0
.
(74)
Therefore, we can directly relate the measured in-
verse decay rate hhΓ(E)
−1
ii
0
to the time-averaged
fluctuations, and from measurement of each for
changing device size or coupling strength, as
shown in the numerical experiments of Section
3, yields information on the scaling of the Hilbert
space dimension. We finally note that it is sim-
ply the lack of direct knowledge of the constant
C
0
β
which prevents measurements where there is
a non-negligible change in the decay rate with
energy (a non-exponential decay to equilibrium)
from constituting a direct measurement of the
value of the Hilbert space dimension. This con-
stant depends on the functional form D(E) and
Γ(E), and thus, if these are known, inference of
the Hilbert space dimension itself is thus obtain-
able.
Finally, we note that the finite temperature ap-
proach above can be extended to the low tem-
perature regime, as shown in Appendix E, from
which we can recover the pure state FDT found
in Ref. [41] in the low temperature limit.
7 Discussion
The results shown above demonstrate how the
chaotic dynamics of thermalization may be ex-
ploited in order to gain information on the com-
plexity of the unitary quantum dynamics of a sys-
tem. We have proposed an experimentally viable
protocol, by which measurements of a local ob-
servable of a probe qubit may be exploited to
measure the Hilbert space dimension of an er-
godic quantum device, initialized in an infinite
temperature state. We note that this measures
the dimension of the states directly involved in
dynamics only, and thus provides a more accurate
measure of the complexity of the dynamics than
a simple estimate of the Hilbert space dimension
from the number of qubits. In this sense, such a
measurement of a large enough quantum device,
if shown to be ergodic in the sense outlined above,
would be a convincing indicator of the so called
‘quantum supremacy’ of the quantum device.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 15
On a practical level, for a generic non-
integrable Hamiltonian, our results may be ob-
served in two ways: measurement of a probe
observable for (i) changing the number of
qubits/ions/... in the quantum device (as in Figs.
3 and 4) , or (ii) changing the probe-bath coupling
(as in Fig. 2). The latter is perhaps the simplest
experimental methodology, which we show can
confirm the ergodic behaviour of a system, that
is, that the unitary dynamics requires an exten-
sive proportion of the Hilbert space, by showing a
linear relationship between the long-time fluctu-
ations and decay rate. For a model where the de-
vice size may be altered, our FDT provides even
deeper insight, allowing also for the experimen-
tal observation of the scaling of the Hilbert space
dimension with system size.
For cases where an exponential decay to equi-
librium is observed, which we show implies that
the decay rate is constant over a large range of
energies, our method allows the experimenter to
access the numerical value of the Hilbert space
dimension itself, not simply its scaling with size
or coupling strength. This can be obtained from
a single time trace of the decay to equilibrium
of the observable, from the measurement of the
decay rate, and fluctuations around equilibrium.
8 Acknowledgements
We acknowledge funding from project PGC2018-
094792-B-I00 (MCIU/AEI/FEDER, UE), EP-
SRC grant no. EP/M508172/1, and from COST
Action CA17113.
A Summary of RMT Approach
A.1 Random Matrix Model
The random matrix model under study may be expressed by a non-interacting part
(H
0
)
αβ
= E
α
δ
αβ
(75)
where E
α
= αω
0
, and ω
0
= 1/N is the spacing between energy levels, and perturbation term, modelled
by a random matrix,
V
αβ
= h
αβ
, (76)
where h
αβ
are independent random numbers selected from the Gaussian Orthogonal Ensemble (GOE),
and rescaled by a coupling strength g, such that the matrix h has the probability distribution,
P (h) ∝ exp
−
N
4g
2
Tr h
2
, (77)
giving hh
αβ
i = 0, and hh
2
αβ
i = g
2
/N for α 6= β, and otherwise hh
2
αα
i = 2g
2
/N .
In Ref. [40] the current authors developed a consistent theoretical model of random wave functions
|ψ
µ
i =
P
α
c
µ
(α)|φ
α
i, for the random matrix model above. We make the ansatz on the probability
distribution on the c
µ
(α)s,
p(c, Λ) =
1
Z
p
e
−
P
µα
c
2
µ
(α)
2Λ(µ,α)
Y
µν
µ>ν
δ(
X
α
c
µ
(α)c
ν
(α)), (78)
where the δ-function term explicitly accounts for the orthogonality of the many-body eigenstates (which
we showed to be necessary in order to obtain a consistent form for the off-diagonal matrix elements).
This distribution Λ(µ, α) can then be shown to be a Lorentzian of width Γ =
πg
2
N ω
0
[40, 45]. From
Eq. (78), one can calculate arbitrary correlation functions of the c
µ
(α) coefficient by first defining the
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 16
generating function,
G
(od)
µν
(
~
ξ
µ
,
~
ξ
ν
) =
Z Z
exp
−
X
α
c
2
µ
(α)
2Λ(µ, α)
+
c
2
ν
(α)
2Λ(ν, α)
+ ξ
µ,α
c
µ
(α) + ξ
ν,α
c
ν
(α)
× δ(
X
α
c
µ
(α)c
ν
(α))
Y
α
dc
µ
(α)dc
ν
(α)
∝ exp
1
2
X
α
ξ
2
µ,α
Λ(µ, α) +
1
2
X
α
ξ
2
ν,α
Λ(ν, α)−
1
2
X
α,β
ξ
µ,α
ξ
µ,β
ξ
ν,α
ξ
ν,β
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
,
(79)
where in the second line we have re-expressed the δ-functions in their Fourier form. The superscript
(od) indicates that this is the ‘off-diagonal’ generating function, requiring µ 6= ν. The diagonal case is
discussed below. The correlation functions may then be calculated by performing successive derivatives
with respect to the force terms ξ via
hc
µ
(α)c
ν
(β) · · · c
µ
(α
0
1
)c
ν
(β
0
1
)i
V
=
1
G
µν
∂
ξ
µ,α
∂
ξ
ν,β
· · · ∂
ξ
µ,α
0
1
∂
ξ
ν,β
0
1
G
µν
ξ
µ,α
=0,ξ
ν,α
=0
.
(80)
In particular, the correlation function hc
µ
(α
0
)c
ν
(β
0
)c
µ
(α)c
ν
(β)i
V
was found in [40] for µ 6= ν to be
equal to
hc
µ
(α
0
)c
ν
(β
0
)c
µ
(α)c
ν
(β)i
V
= Λ(µ, α
0
)Λ(ν, β
0
)δ
α
0
α
δ
β
0
β
−
Λ(µ, α
0
)Λ(ν, α
0
)Λ(µ, α)Λ(ν, α)δ
α
0
β
0
δ
αβ
Λ
(2)
(µ, ν)
−
Λ(µ, α
0
)Λ(ν, α
0
)Λ(µ, β
0
)Λ(ν, β
0
)δ
α
0
β
δ
β
0
α
Λ
(2)
(µ, ν)
,
(81)
with
Λ
(n)
(µ, ν) :=
ω
0
nΓ/π
(E
µ
− E
ν
)
2
+ (nΓ)
2
, (82)
where the superscript (n) is left out for n = 1. The latter two terms in Eq. (81) arise as an explicit
result of the orthogonality factor in Eq. (78).
We stress here that the generating function Eq. (79) explicitly requires µ 6= ν, as it models the
interactions due to mutual orthogonality of two random wave functions. For the diagonal part, we
have the much simpler generating function,
G
(d)
µµ
=
Z
Y
α
dc
µ
(α) exp
"
−
c
2
µ
(α)
2Λ(µ, α)
+ ξ
µ,α
c
µ
(α)
#
∝ exp
"
1
2
X
α
ξ
2
µ,α
Λ(µ, α)
#
(83)
Thus, we have, from successive derivatives with respect to the force term ξ
µ,α
with different non-
interacting indices, and taking the limit as ξ
µ,α
→ 0, as in Eq. (80),
hc
µ
(α)c
µ
(β)c
µ
(α
0
)c
µ
(β
0
)i
V
= Λ(µ, α)Λ(µ, α
0
)δ
αβ
δ
α
0
β
0
+ Λ(µ, α)Λ(µ, β)(δ
αα
0
δ
ββ
0
+ δ
αβ
0
δ
α
0
β
),
(84)
for the diagonal case.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 17
Figure 6: (a) Numerical confirmation of the random matrix FDT for an infinite temperature initial state, Eq. (118)
for observables O
odd
and O
sym
. (b) Shows the random matrix FDT for a high temperature initial state β = 100, and
(c) for a low temperature (β = 500). Here g = 0.04, so Γ ∼ 0.007, and thus the high temperature limit β Γ
−1
is
approximately β 125. We thus observe that the finite temperature result Eqs. (62) and (115) (which we note are
equivalent, the latter is used here), is fulfilled for all temperatures. We note that the low temperature limit above
uses ρ
B
∼ e
−β(E−E
0
)
, with E
0
=
E
max
2
, to ensure that the initial state is not simply the ground state. For this limit
we also use W
O
= [∆O
2
], as discussed in the final section. Simulations are performed with a single realization of
the random matrix V , and thus we observe directly the self-averaging property.
A.2 Self-Averaging
A final property we exploit is that of the self-averaging of random matrices, specifically, in summations
we make the assumption that the observable time dependence and fluctuations are equal to the ensemble
average over perturbations V , such that hO(t)i ≈ hhO(t)ii
V
. This has been previously checked for
this model numerically in Refs. [40, 41], and can also be seen to be valid numerically in this case
from Figs. 6 and 7, which compare our analytical calculations to single realizations of the random
matrix Hamiltonian. Self-averaging is a common tool in RMT [47], and whilst examples of behaviour
that violates a self-averaging assumption indeed exist [48, 49], these are generally due to additional
constraints preventing such statistical behaviour, for global observables such as the survival probability,
or in more exotic regimes e.g. close to critical points [50, 51].
A recent and timely study [48] has looked in detail at the self-averaging behaviour of multiple
quantities for both a full GOE, and a realistic spin chain, finding that in general one should be wary
of simply applying self-averaging. It is shown, for example, that the survival probability is not self-
averaging at any timescale. Indeed, the current authors have recently shown that this quantity cannot
be expected to behave ergodically, and thus RMT should not apply to the survival probability [49].
In the light of such studies, as well as such counterexamples to self-averaging outlined above, it is
important to further justify the use of this property in more detail.
The definition for self-averaging used in [48] is that the vanishing of the quantity
R
O
(t) =
hhO(t)i
2
i
V
− hhO(t)ii
2
V
hhO(t)ii
2
V
, (85)
as the system size increases. Indeed, a very similar quantity was bounded for this RMT model in Ref.
[35]. In Appendix C, we will show that self-averaging in the sense of the vanishing of Eq. (85) indeed
occurs for our model using the assumptions on observables outlined above.
We further note that as the realistic Hamiltonians of interest are not necessarily disordered, in
the sense that they themselves may have no random component, we do not require that they be
self-averaging themselves (indeed this has no real meaning in this case). Rather, the assumption in
our case is that when expressed in that non-interacting basis, a chaotic interaction Hamiltonian is
suitably ‘random’, such that it resembles a single realization of the ensemble we describe above, which
itself is self-averaging. In fact, it is not important that the entire interaction Hamiltonian resembles
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 18
an element of this ensemble, we require instead that locally, for any energy E within the bulk, the
interaction Hamiltonian resembles an element of the random matrix ensemble within a width ∼ Γ from
E.
A.3 RMT Numerics
Here we confirm our analytical results with numerical calculations with the random matrix Hamilto-
nian. In particular, we show in Fig. 6a, that the infinite temperature fluctuation-dissipation theorem
(FDT),
δ
2
O
(∞) =
W
O
ω
0
4πN
B
Γ
,
(86)
is satisfied in this model. This is shown for two ‘observables’ of the RMT model, O
odd
and O
sym
, which
are chosen to be diagonal in the non-interacting basis, with diagonal elements given by,
(O
odd
)
αα
=
(
1 if α = odd
0 otherwise,
(87)
for O
odd
, and
(O
sym
)
αα
=
(
1 if α = odd
−1 otherwise,
(88)
for O
sym
. These observables are chosen, as in Refs. [40, 41], as they resemble realistic observables,
such as local Pauli operators, in the sense that they are well defined, sparse, and highly degenerate
[29] in the non-interacting basis. For our RMT numerical calculations, we define the initial state as
ρ
αβ
= e
−β(E
α
−E
0
)
δ
αβ
δ
α,odd
, (89)
such that hO(0)i = 1. The energy shift E
0
is simply to avoid edge effects at lower temperatures, where
a large fraction of the initial state population would otherwise be in the ground state.
In Fig. 7 we plot the time dependence of the above observables for an infinite temperature initial
state, and compare these to the observable time dependence of Eq. (37), derived below.
The high temperature limit, in which our FDT is derived, is defined by β
−1
Γ. In the numerics,
we use g = 0.05, and thus Γ ∼ 0.007, so the high temperature limit requires β . 125. We show plots
for β = 100 and β = 500 in Figs. 6b and 6c, respectively. For the parameters used these correspond
to a high temperature, near the edge of the expected limit, and a low temperature initial state. We
observe that the finite temperature form in fact works well for all β values. This is further discussed
analytically in the final section below, where we see that the low temperature limit of our current
approach is equal to the pure state result previously obtained in Ref. [41].
B Gaussian, Non-Gaussian, and Mixed Contractions
Here we calculate the contributions to the long-time fluctuations of Eq. (47). We see that the relevant
correlation function is an 8-point correlation function, which we will show below may be split into
three groups: products of Gaussian contractions, products of non-Gaussian contractions, and mixed
products of two Gaussian and one non-Gaussian contraction.
An example of the first form, Gaussian contractions only, is
hc
µ
(α)c
ν
(α)c
µ
(β)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
)c
µ
(β
0
)c
ν
(β
0
)i
V
= Λ(µ, α)Λ(ν, α)Λ(µ, α
0
)Λ(ν, α
0
)δ
αβ
δ
α
0
β
0
. (90)
Indeed, there are three such contractions, occurring each time between pairs of indices, that only
contribute two kronecker-δ factors - ∼ δ
αβ
δ
α
0
β
0
, δ
αα
0
δ
ββ
0
, δ
αβ
0
δ
βα
0
. Other Gaussian contractions may
be defined, but may be ignored due to a reduction in the number of summations.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 19
Figure 7: Time dependence of random matrix observables O
odd
and O
sym
. Exact diagonalization numerics (solid
lines) show time evolutions for a single realization of the random matrix perturbation V . RMT calculation (dashed
lines), show Eq. (37), with hO(t)i
0
= 1, and hOi
MC
= [O
αα
] = 0(0.5) for O = O
sym
(O
odd
). Parameters used are
N = 500, g = 0.05, β = 0, 100.
In a similar manner, we may define non-Gaussian contractions that do not reduce the number of
summations at all, such that they contribute on the same order as the Gaussian contractions above.
An example is,
hc
µ
(α)c
ν
(α)c
µ
(β)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
)c
µ
(β
0
)c
ν
(β
0
)i
V
= L
ααββ
µν
L
α
0
α
0
β
0
β
0
µν
. (91)
It can be easily seen that there are three non-Gaussian contractions of this form, with the other two
being defined by swapping pairs of primed and unprimed indices in turn.
Finally, we see that mixed contractions may also contribute, for example,
hc
µ
(α)c
ν
(α)c
µ
(β)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
)c
µ
(β
0
)c
ν
(β
0
)i
V
= Λ(µ, α)Λ(ν, α)L
α
0
α
0
β
0
β
0
µν
δ
αβ
. (92)
We can see that there are six terms of this form that contribute only one δ factor. These are
δ
αβ
, δ
αα
0
, δ
αβ
0
, δ
βα
0
, δ
ββ
0
, δ
α
0
β
0
.
We thus obtain that for the contribution from Gaussian contractions, δ
2
G
(∞),
δ
2
G
(∞) =
X
µν
µ6=ν
X
αβα
0
β
0
w
α
w
β
O
α
0
α
0
O
β
0
β
0
Λ(µ, α)Λ(ν, β)Λ(µ, α
0
)Λ(ν, β
0
)(δ
αβ
δ
α
0
β
0
+ δ
αα
0
δ
ββ
0
+ δ
αβ
0
δ
α
0
β
)
=
X
µν
µ6=ν
X
αβ
Λ(µ, α)Λ(ν, α)Λ(µ, β)Λ(ν, β)[w
2
α
O
2
ββ
+ 2w
α
w
β
O
αα
O
ββ
]
(93)
Similarly, for the non-Gaussian (δ
2
NG
(∞)), and mixed (δ
2
M
(∞)) contractions we have
δ
2
NG
(∞) = 3
X
µν
µ6=ν
X
αβα
0
β
0
w
α
w
β
O
α
0
α
0
O
β
0
β
0
L
αβαβ
µν
L
α
0
β
0
α
0
β
0
µν
,
(94)
and
δ
2
M
(∞) = −
X
µν
µ6=ν
X
αα
0
β
0
Λ(µ, α)Λ(ν, α)L
α
0
β
0
α
0
β
0
µν
w
2
α
O
α
0
α
0
O
β
0
β
0
+ 4w
α
O
αα
w
α
0
O
β
0
β
0
+ O
2
αα
w
α
0
w
β
0
,
(95)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 20
respectively. In order to perform the summations over non-interacting indices in Eq. (47) we define
course grained averages of observable elements O
αα
as in Eq. (23). The key assumption in writing
(23) is that the average,
[O
αα
]
µ
:=
X
α
Λ(µ, α)O
αα
, (96)
changes slowly in energy E
µ
with respect to the width Γ. Similar averages over the initial state, or
mixed averages must also be defined, such as
X
α
Λ(µ, α)Λ(ν, α)w
α
= [w
α
]
µ
Λ
(2)
(µ, ν)
X
α
Λ(µ, α)Λ(ν, α)w
α
O
αα
= [w
α
O
αα
]
µ
Λ
(2)
(µ, ν).
(97)
Note that, [O
αα
]
µ
can be interpreted as a microcanonical average of the observable O.
Now, using Eqs. (23) and (97), we obtain,
δ
2
G
(∞) =
X
µν
µ6=ν
[w
2
α
]
µ
[O
2
αα
]
µ
+ 2[w
α
O
αα
]
2
µ
Λ
(2)
(µ, ν)
2
,
δ
2
NG
(∞) = 3
X
µν
µ6=ν
[w
α
]
2
µ
[O
αα
]
2
µ
Λ
(2)
(µ, ν)
2
,
δ
2
M
(∞) = −
X
µν
µ6=ν
Λ
(2)
(µ, ν)
2
[w
2
α
]
µ
[O
αα
]
2
µ
+ 4[w
α
]
µ
[O
αα
]
µ
[w
α
O
αα
]
µ
+ [w
α
]
2
µ
[O
2
αα
]
µ
.
(98)
C Proof of Self-Averaging
In addition to the numerical demonstration in Figs. 6 and 7, we are now able to analytically demon-
strate that self-averaging, in the sense of Eq. (85), occurs for this model. We have already obtained
hhO(t)ii
V
in Section 5 (note that here we assumed self-averaging, and thus simply used the notation
hO(t)i, in the current section we require the self-averaging step to be explicit, so use the angle brackets
h· · · i
V
whenever self-averaging is performed), and thus, what remains is a calculation of hhO(t)i
2
i
V
.
We can do this by writing,
∆O(t)
2
=
X
µν
µ6=ν
X
µ
0
ν
0
µ
0
6=ν
0
X
αβα
0
β
0
X
α
1
β
1
α
0
1
β
0
1
w
αβ
w
α
0
β
0
O
α
0
β
0
O
α
0
1
β
0
1
× c
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(β
0
)c
µ
0
(α
1
)c
ν
0
(β
1
)c
µ
0
(α
0
1
)c
ν
0
(β
0
1
)
=
X
µν
µ6=ν
X
µ
0
ν
0
µ
0
6=ν
0
X
αβα
0
X
α
1
β
1
α
0
1
N
O
X
n,n
1
w
αβ
w
α
0
β
0
O
α
0
,α
0
+n
O
α
0
1
,α
0
1
+n
1
× c
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
+ n)c
µ
0
(α
1
)c
ν
0
(β
1
)c
µ
0
(α
0
1
)c
ν
0
(α
0
1
+ n
1
)
(99)
from which, we see that the relevant correlation function after self-averaging is hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
+
n)c
µ
0
(α
1
)c
ν
0
(β
1
)c
µ
0
(α
0
1
)c
ν
0
(α
0
1
+ n
1
)i
V
. Indeed, the contribution from Gaussian contractions can easily
be seen to be identical to that of h∆O(t)i
2
V
, as these correlators only contract identical interacting
indices µ, ν. Furthermore, we see that when non-Gaussian contractions are defined between primed
and unprimed interacting indices, µ, ν, there is a reduction in the number of summations. For example,
we compare the contribution
X
µνµ
0
ν
0
X
αβα
0
X
α
1
β
1
α
0
1
hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
)i
V
hc
µ
0
(α
1
)c
ν
0
(β
1
)c
µ
0
(α
0
1
)c
ν
0
(α
0
1
)i
V
, (100)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 21
to that with mixed interacting indices,
X
µνµ
0
ν
0
X
αβα
0
X
α
1
β
1
α
0
1
hc
µ
(α)c
ν
0
(β
1
)c
µ
(α
0
)c
ν
0
(α
0
1
)i
V
hc
µ
0
(α
1
)c
ν
(β)c
µ
0
(α
0
1
)c
ν
(α
0
)i
V
, (101)
and observe that the latter contribution is negligible, due to the extra reduced summation, as there are
no repeated non-interacting indices within the non-Gaussian contraction. We note that the dominating
contribution from non-Gaussian contractions is inevitably from the n, n
1
= 0 term shown above (as
this reduced summation does not change the order of the contribution as N
O
N ). Thus we have
a single dominating contribution given by the 4-point correlation functions defined with primed, or
unprimed, interaction indices only:
h∆O(t)
2
i
V
=
X
µν
µ6=ν
X
µ
0
ν
0
µ
0
6=ν
0
X
αβα
0
X
α
1
β
1
α
0
1
N
O
X
n,n
1
w
αβ
w
α
0
β
0
O
α
0
,α
0
+n
O
α
0
1
,α
0
1
+n
1
e
−i(E
µ
−E
ν
+E
µ
0
−E
ν
0
)t
× hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
+ n)c
µ
0
(α
1
)c
ν
0
(β
1
)c
µ
0
(α
0
1
)c
ν
0
(α
0
1
+ n
1
)i
V
=
X
µν
µ6=ν
X
αβα
0
N
O
X
n
w
αβ
O
α
0
,α
0
+n
e
−i(E
µ
−E
ν
)t
hc
µ
(α)c
ν
(β)c
µ
(α
0
)c
ν
(α
0
+ n)i
V
×
X
µ
0
ν
0
µ
0
6=ν
0
X
α
1
β
1
α
0
1
N
O
X
n
1
w
α
0
β
0
O
α
0
1
,α
0
1
+n
1
e
−i(E
µ
0
−E
ν
0
)t
hc
µ
0
(α
1
)c
ν
0
(β
1
)c
µ
0
(α
0
1
)c
ν
0
(α
0
1
+ n
1
)i
V
= h∆O(t)i
2
V
.
(102)
Therefore, we can see that the transient component ∆O(t) of the time evolution is indeed self-averaging.
Now, we can similarly see that the long-time average is self-averaging via,
O
MC
=
X
µ
X
αβ
w
αβ
O
µµ
c
µ
(α)c
µ
(β)
=
X
µ
X
αβα
0
β
0
w
αβ
O
α
0
β
0
c
µ
(α)c
µ
(β)c
µ
(α
0
)c
µ
(β
0
)
=
X
µ
X
αβα
0
N
O
X
n
w
αβ
O
α
0
,α
0
+n
c
µ
(α)c
µ
(β)c
µ
(α
0
)c
µ
(α
0
+ n).
(103)
The ensemble average of this is then,
hO
MC
i
V
=
X
µ
X
αβα
0
N
O
X
n
w
αβ
O
α
0
,α
0
+n
hc
µ
(α)c
µ
(β)c
µ
(α
0
)c
µ
(α
0
+ n)i
V
=
X
µ
X
αα
0
w
αα
O
α
0
α
0
Λ(µ, α)Λ(µ, α
0
),
(104)
where one can easily see that the terms with n 6= 0 require an additional contraction, and thus the
contribution is on a lower order, with fewer summations over the non-interacting indices. We then see
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 22
that,
O
2
MC
=
X
µν
X
αβα
1
β
1
w
αβ
w
α
1
β
1
O
µµ
O
νν
c
µ
(α)c
µ
(β)c
ν
(α
1
)c
ν
(β
1
)
=
X
µν
X
αβα
1
β
1
X
α
0
β
0
α
0
1
β
0
1
w
αβ
w
α
1
β
1
O
α
0
β
0
O
α
0
1
β
0
1
c
µ
(α)c
µ
(β)c
ν
(α
1
)c
ν
(β
1
)c
µ
(α
0
)c
µ
(β
0
)c
ν
(α
0
1
)c
ν
(β
0
1
)
=
X
µν
X
αβα
1
β
1
X
α
0
α
0
1
N
O
X
nn
0
w
αβ
w
α
1
β
1
O
α
0
,α
0
+n
O
α
0
1
,α
0
1
+n
0
× c
µ
(α)c
µ
(β)c
ν
(α
1
)c
ν
(β
1
)c
µ
(α
0
)c
µ
(α
0
+ n)c
ν
(α
0
1
)c
ν
(α
0
1
+ n
0
),
(105)
which, after self-averaging becomes,
hO
2
MC
i
V
=
X
µν
X
αβα
1
β
1
X
α
0
α
0
1
N
O
X
nn
0
w
αβ
w
α
1
β
1
O
α
0
,α
0
+n
O
α
0
1
,α
0
1
+n
0
× hc
µ
(α)c
µ
(β)c
ν
(α
1
)c
ν
(β
1
)c
µ
(α
0
)c
µ
(α
0
+ n)c
ν
(α
0
1
)c
ν
(α
0
1
+ n
0
)i
V
=
X
µν
X
αα
1
X
α
0
α
0
1
w
αα
w
α
1
α
1
O
α
0
α
0
O
α
0
1
α
0
1
Λ(µ, α)Λ(ν, α
1
)Λ(µ, α
0
)Λ(ν, α
0
1
)
= hO
MC
i
2
V
(106)
where we see once more that the n, n
0
= 0 terms dominate, and that contractions between the primed
and unprimed indices induce additional kronecker-delta factors, and thus are negligible.
Thus, we have that R = 0 (see Eq. (85)), and we analytically observe self-averaging in this model.
We note here that in the derivation of the 4-point correlator, we assume that the only ‘interactions’
between the eigenstates themselves are those enforced due to the mutual orthogonality of eigenstates
via the condition
P
α
c
µ
(α)c
ν
(α) = δ
µν
[40]. As these interactions occur pairwise between eigenstates,
correlation functions with more than two interacting indices µ, ν contribute an additional δ
µν
.
D Bound of dynamical term
In this section we bound third term in Eq. (33), showing that it’s contribution is negligible. To do so,
we proceed by defining,
A(t) =
X
µν
µ6=ν
X
αβ
w
αβ
O
αβ
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
e
−i(E
µ
−E
ν
)t
. (107)
We may now use the relation |
P
i
a
i
| ≤
P
i
|a
i
|, to write
|A(t)| ≤
X
µν
µ6=ν
X
αβ
w
αβ
O
αβ
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
e
−i(E
µ
−E
ν
)t
≤
X
µν
µ6=ν
X
αβ
w
αβ
O
αβ
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
=
X
αβ
|w
αβ
O
αβ
|
X
µν
µ6=ν
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
≤
3ω
0
4πΓ
X
α
|w
αβ
O
αβ
| ,
(108)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 23
where we have used that,
X
µν
µ6=ν
Λ(µ, α)Λ(µ, β)Λ(ν, α)Λ(ν, β)
Λ
(2)
(µ, ν)
= ω
0
(E
α
− E
β
)
2
Γ + 12Γ
3
π((E
α
− E
β
)
2
+ 4Γ
2
)
2
≤
3ω
0
4πΓ
. (109)
Now, using
P
αβ
O
αβ
=
P
α
P
N
O
n
O
α,α+n
δ
β,α+n
[41], we have
X
αβ
|w
αβ
O
αβ
| =
X
α,n
|w
α,α+n
O
αα+n
|
≤ max
αβ
(O
αβ
)
X
α,n
w
α,α+n
.
(110)
Now, we see that in the case of w
αβ
∼ δ
αβ
, the bound simply becomes
|A(t)| ≤ max
α
(O
αα
)
3ω
0
4πΓ
, (111)
similarly, this is the case for our application studied in the main text, where O
αβ
∼ δ
αβ
. We note that
the condition of diagonal w
αβ
correspond to a reasonable initial state in many experimental situations,
since thermalization takes place typically by incoherent exchange of energy in the basis of eigenstates
of H
0
. For states with coherences, to bound this quantity we require that the coherences are not large
on the off-diagonals defined by α, α + n, in the sense that
P
α
w
α,α+n
. O(1), so
|A(t)| . N
O
max
α
(O
αα
)
3ω
0
4πΓ
. (112)
In this sense, ‘special’ initial states may be chosen that do not satisfy this bound, but they are highly
atypical. Further, we note that the stronger of the two assumptions made on the form of observables,
that is that we may write
P
αβ
O
αβ
=
P
α
P
N
O
n
O
α,α+n
δ
β,α+n
, is only used above in bounding |A(t)|,
the other terms are more general. As noted above, whilst this form is reasonable for generic local
variables [41], one may of course be able to build observables that do not match this form, in which
case, one may either have an additional contribution due to A(t), or be able to arrive at a similar
bound.
E FDT in terms of Thermal Averages
In this section, we show that the FDT may be recast in terms of the thermal averages hh· · · ii
β
, obtained
from the a fit to the decay, shown in Sec 5. We begin by re-expressing our finite temperature FDT in
terms of the thermal averages,
hhA(E)ii
β
:=
1
Z
0
β
Z
∆E
0
dED(E)e
−β(E−E
0
)
A(E), (113)
with Z
0
β
:=
R
∆E
0
dED(E)e
−β(E−E
0
)
, where we have introduced the low-energy cut-off E
0
, which is the
energy of our zero temperature pure state. The addition of this cut-off ensures that, for non-zero E
0
,
the initial state is not the ground state of the bath at zero temperature.
To include this thermal average, we return to Eq. (70), and note that,
Z
∆E
0
dE
e
−2βE
Γ(E)
= Z
0
2β
hhD(E)
−1
Γ(E)
−1
ii
2β
. (114)
We thus obtain,
δ
2
O
(∞) =
W
O
Z
0
2β
8πZ
2
β
hhD(E)
−1
Γ(E)
−1
ii
2β
.
(115)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 24
We note that this can be put in the same form as the infinite temperature case, δ
2
O
(∞) ∼ Γ
−1
, by
noting that the quantity
C
0
β
=
hhD(E)
−1
Γ(E)
−1
ii
2β
hhD(E)ii
−1
2β
hhΓ(E)
−1
ii
2β
(116)
depends only on the particular forms of the functions D(E) and Γ(E), and the temperature - impor-
tantly, not on N, or on the system-bath coupling strengths (for weak couplings), as both the numerator
and denominator share the same dependence in N and the system-bath coupling in the thermodynamic
limit. As such, we can write
δ
2
O
(∞) = C
0
β
W
O
Z
0
2β
8πZ
2
β
hhD(E)ii
2β
hhΓ(E)
−1
ii
2β
,
(117)
and thus we recover the form of our main result, Eq. (3) of the main text, with χ(N) =
C
0
β
W
O
Z
2β
16πZ
2
β
hhD(E)ii
β
.
For the random matrix case, where D(E)
−1
= ω
0
, and Γ(E) = Γ is constant, we also have,
δ
2
O
(∞) =
W
O
Z
0
2β
ω
0
8πZ
2
β
Γ
(118)
as in this case, the thermal average hh· · · ii
β
and unbiased thermal average h· · · i
β
, can be seen to be
equal.
We can check that, as required, one obtains the infinite temperature limit derived above by sending
β → 0 by noting that for the infinite temperature case, we have
δ
2
O
(∞) =
W
O
2N
2
B
Z
∆E
0
dE
1
4πΓ(E)
, (119)
which, in terms of the infinite temperature thermal average, may be written (noting that lim
β→0
Z
0
β
=
2N
B
),
δ
2
O
(∞) =
W
O
4πN
B
hhD(E)
−1
Γ(E)
−1
ii
0
= C
0
0
W
O
4πN
B
hhD(E)ii
0
hhΓ(E)
−1
ii
0
.
(120)
We finally note here that it is the measurement of a thermal average of the decay rate, and not the
unbiased average, that obscures the direct measurement of N
B
, such that in general a FDT of the
form of Eq. (120) would be observed experimentally. We thus limit our claim in cases where the
energy dependence of the decay rate is important (which, as can be seen from Eq. (40), is implied by
a non-exponential decay to equilibrium, due to multiple contributing decay rates) to the experimental
verification of the finite size scaling of the Hilbert space dimension, rather than the particular value of
N
B
, without additional assumptions or the numerical calculation of C
0
β
.
E.1 Low Temperature FDT
We now turn to the low temperature limit of Eq. (67) for which we expect to obtain the same result
as the pure state case of Ref. [41], given by,
δ
2
O
(∞) =
[∆O
2
αα
]
4πD(E
α
0
)Γ
, (121)
where D(E
α
0
) is the DOS at the initial state energy E
α
0
, which is chosen to be in the bulk of the
spectrum, and [∆O
2
αα
] := [O
2
αα
] − [O
αα
]
2
. We have that in this limit,
hhA(E)ii
∞
= A(E
0
), (122)
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 25
so
C
0
∞
=
D(E
0
)
−1
Γ(E
0
)
−1
D(E
0
)
−1
Γ(E
0
)
−1
= 1, (123)
and thus,
δ
2
O
(∞) =
W
O
4πD(E
0
)
Γ(E
0
)
−1
,
(124)
which is the zero temperature limit, Eq. (121) when W
0
= [∆O
2
αα
]. Which can be seen to be the case
for zero temperature as follows. Recalling that W
0
is defined by
W
O
=
W
µ
[w
2
α
]
(125)
where,
W
µ
= [w
2
α
]
µ
[O
2
αα
]
µ
+ 2[w
α
O
αα
]
2
µ
+ 3[w
α
]
2
µ
[O
αα
]
2
µ
− [w
2
α
]
µ
[O
αα
]
2
µ
− 4[w
α
]
µ
[O
αα
]
µ
[w
α
O
αα
]
µ
− [w
α
]
2
µ
[O
2
αα
]
µ
.
(126)
We see that in the zero temperature limit w
α
∼ δ
αα
0
, and thus, the averages in Eq. (97) contribute
to a lower order as the number of summations is reduced for terms with, e.g. w
α
w
β
, over terms with,
e.g. w
2
α
. This can be seen to lead to [w
α
]
2
µ
[w
2
α
]
µ
. Similarly, both terms above with mixed averages
[w
α
O
αα
]
µ
contribute on the order [w
α
]
2
µ
, as [w
α
O
αα
]
µ
= [w
α
]
µ
O
↑
. Using only the remaining terms, we
have that,
W
µ
= [w
2
α
]
µ
[O
2
αα
]
µ
− [w
2
α
]
µ
[O
αα
]
2
µ
= [w
2
α
]
µ
[∆O
2
αα
],
(127)
so,
W
O
= [∆O
2
], (128)
as required.
We recall that until Eq. (67), no assumptions on the initial state or observable are made other than
the ability to define the required microcanonical averages. As such, taking the low temperature limit
at this point, as we have done above, does not contradict any assumptions made.
References
[1] Immanuel Bloch, Jean Dalibard, and Sylvain Nascimbène. Quantum simulations with ultracold
quantum gases. 8(4):267–276, 2012. DOI: https://doi.org/10.1038/nphys2259.
[2] Maciej Lewenstein, Anna Sanpera, and Veronica Ahufinger. Ultracold Atoms in Optical
Lattices: Simulating quantum many-body systems. Oxford University Press, 2012. DOI:
https://doi.org/10.1080/00107514.2013.800135.
[3] M Aidelsburger, J. L. Ville, R Saint-Jalm, S. Nascimbène, J Dalibard, and J Beugnon. Relaxation
Dynamics in the Merging of N Independent Condensates. Phys. Rev. Lett., 119(19):190403, 2017.
DOI: https://doi.org/10.1103/PhysRevLett.119.190403.
[4] D. Porras and J. I. Cirac. Effective quantum spin systems with trapped ions. Phys. Rev. Lett.,
92:207901, May 2004. DOI: https://doi.org/10.1103/PhysRevLett.92.207901.
[5] Ch Schneider, Diego Porras, and Tobias Schaetz. Experimental quantum simulations of
many-body physics with trapped ions. Rep. Prog. Phys., 75(2):024401, jan 2012. DOI:
https://doi.org/10.1088/0034-4885/75/2/024401.
[6] R Blatt and C F Roos. Quantum simulations with trapped ions. Nature Physics, 8(4):277–284,
2012. DOI: https://doi.org/10.1038/nphys2252.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 26
[7] Govinda Clos, Diego Porras, Ulrich Warring, and Tobias Schaetz. Time-Resolved Observation
of Thermalization in an Isolated Quantum System. Phys. Rev. Lett., 117(17):1–6, 2016. DOI:
https://doi.org/10.1103/PhysRevLett.117.170401.
[8] Hyosub Kim, Yeje Park, Kyungtae Kim, H. S. Sim, and Jaewook Ahn. Detailed Balance of
Thermalization dynamics in Rydberg atom quantum simulators. Phys. Rev. Lett., 120(18):180502,
2018. DOI: https://doi.org/10.1103/PhysRevLett.120.180502.
[9] Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes
Pichler, Soonwon Choi, Alexander S Zibrov, Manuel Endres, Markus Greiner, Vladan Vuletic,
and Mikhail D Lukin. Probing many-body dynamics on a 51-atom quantum simulator. Nature,
551(7682):579–584, 2017. DOI: https://doi.org/10.1038/nature24622.
[10] C Neill, P Roushan, M Fang, Y. Chen, M Kolodrubetz, Z Chen, A Megrant, R Barends, B Camp-
bell, B Chiaro, A Dunsworth, E. Jeffrey, J Kelly, J Mutus, P. J.J. O’Malley, C Quintana, D Sank,
A Vainsencher, J Wenner, T C White, A Polkovnikov, and J M Martinis. Ergodic dynamics and
thermalization in an isolated quantum system. Nature Physics, 12(11):1037–1041, 2016. DOI:
https://doi.org/10.1038/nphys3830.
[11] P Roushan, E Lucero, John M Martinis, B Chiaro, A Megrant, K Kechedzhi, A Dunsworth,
J Wenner, P Klimov, B Burkett, K Arya, A Vainsencher, J Mutus, H Neven, A Fowler, Z Chen,
Y. Chen, R Barends, S V Isakov, M Giustina, T Huang, J Kelly, M Neeley, T C White, S Boixo,
D Sank, B Foxen, V Smelyanskiy, R Graff, E Jeffrey, C Quintana, and C Neill. A blueprint
for demonstrating quantum supremacy with superconducting qubits. Science, 360(6385):195–199,
2018. DOI: https://doi.org/10.1126/science.aao4309.
[12] Artur S L Malabarba, Luis Pedro Garcia-Pintos, Noah Linden, Terence C. Farrelly, and Anthony J
Short. Quantum systems equilibrate rapidly for most observables. Phys. Rev. E, 90(1), 2014. DOI:
https://doi.org/10.1103/PhysRevE.90.012121.
[13] Luis Pedro García-Pintos, Noah Linden, Artur S.L. Malabarba, Anthony J Short, and Andreas
Winter. Equilibration time scales of physically relevant observables. Phys. Rev. X, 7(3), 2017.
DOI: https://doi.org/10.1103/PhysRevX.7.031027.
[14] Jonas Richter, Jochen Gemmer, and Robin Steinigeweg. Impact of eigenstate ther-
malization on the route to equilibrium. Phys. Rev. E, 99:050104, May 2019. DOI:
https://doi.org/10.1103/PhysRevE.99.050104.
[15] Thiago R. De Oliveira, Christos Charalambous, Daniel Jonathan, Maciej Lewenstein, and Arnau
Riera. Equilibration time scales in closed many-body quantum systems. New J. Phys., 20(3):
33032, 2018. DOI: https://doi.org/10.1088/1367-2630/aab03b.
[16] Fausto Borgonovi, Felix M Izrailev, and Lea F Santos. Exponentially fast dynamics in the
Fock space of chaotic many-body systems. Phys. Rev. E, 99(1):010101(R), 2019. DOI:
https://doi.org/10.1103/PhysRevE.99.010101.
[17] Mauro Schiulaz, E Jonathan Torres-Herrera, and Lea F Santos. Thouless and Relaxation
Time Scales in Many-Body Quantum Systems. Phys. Rev. B, 99(17):174313, 2019. DOI:
https://doi.org/10.1103/PhysRevB.99.174313.
[18] Álvaro M. Alhambra, Jonathon Riddell, and Luis Pedro García-Pintos. Time evolution of correla-
tion functions in quantum many-body systems. 2019. URL http://arxiv.org/abs/1906.11280.
[19] Mark Srednicki. Chaos and quantum thermalization. Phys. Rev. E, 50(2):888–901, 1994. DOI:
https://doi.org/10.1103/PhysRevE.50.888.
[20] Mark Srednicki. The approach to thermal equilibrium in quantized chaotic systems. J. Phys. A:
Math. Gen., 32(3299):1163–1175, 1999. DOI: https://doi.org/10.1088/0305-4470/32/7/007.
[21] Marcos Rigol, Vanja Dunjko, and Maxim Olshanii. Thermalization and its mecha-
nism for generic isolated quantum systems. Nature, 452(7189):854–858, 2008. DOI:
https://doi.org/10.1038/nature06838.
[22] Peter Reimann. Foundation of statistical mechanics under experimentally realistic conditions.
Phys. Rev. Lett., 101(19), 2008. DOI: https://doi.org/10.1103/PhysRevLett.101.190403.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 27
[23] V. I. Yukalov. Equilibration and thermalization in finite quantum systems. Laser Physics Letters,
8(7):485–507, 2011. DOI: https://doi.org/10.1002/lapl.201110002.
[24] Marcos Rigol and Mark Srednicki. Alternatives to eigenstate thermalization. Phys. Rev. Lett.,
108(11):1–5, 2012. DOI: https://doi.org/10.1103/PhysRevLett.108.110601.
[25] 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
(3):239–362, 2016. DOI: https://doi.org/10.1080/00018732.2016.1198134.
[26] Wouter Beugeling, R Moessner, and Masudul Haque. Finite-size scaling of eigenstate thermaliza-
tion. Phys. Rev. E, 89(4), 2014. DOI: https://doi.org/10.1103/PhysRevE.89.042112.
[27] Wouter Beugeling, Roderich Moessner, and Masudul Haque. Off-diagonal matrix elements
of local operators in many-body quantum systems. Phys. Rev. E, 91(1), 2015. DOI:
https://doi.org/10.1103/PhysRevE.91.012144.
[28] J Eisert, M Friesdorf, and C Gogolin. Quantum many-body systems out of equilibrium. Nature
Physics, 11(2):124–130, 2014. DOI: https://doi.org/10.1038/nphys3215.
[29] Fabio Anza, Christian Gogolin, and Marcus Huber. Eigenstate Thermalization for Degenerate Ob-
servables. Phys. Rev. Lett., 120(15), 2018. DOI: https://doi.org/10.1103/PhysRevLett.120.150603.
[30] Joshua M Deutsch. Eigenstate thermalization hypothesis. Rep. Prog. Phys., 81(8):082001, 2018.
DOI: https://doi.org/10.1088/1361-6633/aac9f1.
[31] Zeeya Merali. The new thermodynamics: how quantum physics is bending the rules. Nature, 551:
20–22, 2017. DOI: https://doi.org/10.1038/551020a.
[32] Marcos Rigol. Quantum quenches and thermalization in one-dimensional fermionic systems. Phys.
Rev. A, 80(5), 2009. DOI: https://doi.org/10.1103/PhysRevA.80.053607.
[33] Peter Reimann. Typical fast thermalization processes in closed many-body systems. Nature
Communications, 7, 2016. DOI: https://doi.org/10.1038/ncomms10821.
[34] Fausto Borgonovi, F. M. Izrailev, Lea F Santos, and V. G. Zelevinsky. Quantum chaos and
thermalization in isolated systems of interacting particles. Physics Reports, 626:1–58, 2016. DOI:
https://doi.org/10.1016/j.physrep.2016.02.005.
[35] Lennart Dabelow and Peter Reimann. Perturbed relaxation of quantum many-body systems.
2019. URL http://arxiv.org/abs/1903.11881.
[36] Peter Reimann and Lennart Dabelow. Typicality of Prethermalization. Phys. Rev. Lett., 122(8),
2019. DOI: https://doi.org/10.1103/PhysRevLett.122.080603.
[37] Eduardo Jonathan Torres-Herrera, Jonathan Karp, Marco Tavora, and Lea F Santos. Realistic
many-body quantum systems vs. full random matrices: Static and dynamical properties. Entropy,
18(10), 2016. DOI: https://doi.org/10.3390/e18100359.
[38] Ryusuke Hamazaki and Masahito Ueda. Random-matrix behavior of quantum nonintegrable
many-body systems with dyson’s three symmetries. Phys. Rev. E, 99:042116, 2019. DOI:
https://doi.org/10.1103/PhysRevE.99.042116.
[39] R Kubo. The fluctuation-dissipation theorem. Rep. Prog. Phys., 29(1):255, 1966. DOI:
https://doi.org/10.1088/0034-4885/29/1/306.
[40] Charlie Nation and Diego Porras. Off-diagonal observable elements from random matrix theory:
Distributions, fluctuations, and eigenstate thermalization. New J. Phys., 20(10):103003, 2018.
DOI: https://doi.org/10.1088/1367-2630/aae28f.
[41] Charlie Nation and Diego Porras. Quantum chaotic fluctuation-dissipation theorem: Effec-
tive Brownian motion in closed quantum systems. Phys. Rev. E, 99(5):052139, 2019. DOI:
https://doi.org/10.1103/PhysRevE.99.052139.
[42] Phillip Weinberg and Marin Bukov. QuSpin: a Python Package for Dynamics and Exact Diago-
nalisation of Quantum Many Body Systems part I: spin chains. SciPost Phys, 2:003, 2016. DOI:
https://doi.org/10.21468/SciPostPhys.2.1.003.
[43] Phillip Weinberg and Marin Bukov. QuSpin: a Python Package for Dynamics and Exact Diago-
nalisation of Quantum Many Body Systems part I: spin chains. SciPost Phys., 2:003, 2017. DOI:
https://doi.org/10.21468/SciPostPhys.2.1.003.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 28
[44] J. M. Deutsch. Quantum statistical mechanics in a closed system. Phys. Rev. A, 43(4):2046–2049,
1991. DOI: https://doi.org/10.1103/PhysRevA.43.2046.
[45] J. M. Deutsch. A closed quantum system giving ergodicity. (unpublished), 1991. URL https:
//deutsch.physics.ucsc.edu/pdf/quantumstat.pdf.
[46] Peter Reimann. Eigenstate thermalization: Deutsch’s approach and beyond. New J. Phys., 17(5),
2015. DOI: https://doi.org/10.1088/1367-2630/17/5/055025.
[47] Ralf R Müller. Random matrices, free probability and the replica method. In European Signal
Processing Conference, volume 06-10-Sept, pages 189–196, 2015. ISBN 9783200001657. URL
https://ieeexplore.ieee.org/document/7079688.
[48] Mauro Schiulaz, E Jonathan Torres-Herrera, Francisco Pérez-Bernal, and Lea F Santos. Self-
averaging in many-body quantum systems out of equilibrium. 2019. URL http://arxiv.org/
abs/1906.11856.
[49] Charlie Nation and Diego Porras. Non-ergodic quantum thermalization. 2019. URL https:
//arxiv.org/pdf/1908.11773.pdf.
[50] Amnon Aharony and A Brooks Harris. Absence of self-averaging and universal fluctuations
in random systems near critical points. Phys. Rev. Lett., 77(18):3700–3703, 1996. DOI:
https://doi.org/10.1103/PhysRevLett.77.3700.
[51] Giorgio Parisi and Nicolas Sourlas. Scale Invariance in Disordered Systems: The Ex-
ample of the Random-Field Ising Model. Phys. Rev. Lett., 89(25), 2002. DOI:
https://doi.org/10.1103/PhysRevLett.89.257204.
Accepted in Quantum 2019-10-29, click title to verify. Published under CC-BY 4.0. 29