Ergodicity probes: using time-ﬂuctuations 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 ﬂuctuation-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-ﬂuctuations.
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 [13], trapped ions
Charlie Nation: C.Nation@sussex.ac.uk
Diego Porras: D.Porras@iﬀ.csic.es
[47], Rydberg atoms [8, 9], and superconducting
circuits[10, 11]. These systems can be used for
the quantum simulation of many-body models,
or diﬀerent 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 eﬃcient way
can be useful in the diagnosis and certiﬁcation 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 eﬀect of disor-
der or the possible lack of connectivity between
diﬀerent 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 eﬀective
Hilbert space dimension that more accurately de-
scribes the complexity of the system, in terms of
some eﬀective 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 [1218] 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 eﬀective number of states that contribute
to the dynamics of a local observable. Indeed,
above have inspired a bounty of theoretical work
in the ﬁeld of quantum thermalization [1930]. 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
ﬁeld [31].
We assume a quantum quench scenario [8, 21,
28, 32] in which a quantum system is initialized
in a fully-decohered, inﬁnite 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, 3338]. We show that the time-
ﬂuctuations 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
inﬁnite-temperature ﬂuctuation dissipation the-
orem (FDT) [39] for the dynamics of the probe
qubit in order to extract information on the scal-
ing of the eﬀective 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 ﬁnite temperatures, in Section 6.
Finally, we summarise our ﬁndings in Section 7.
Further details and proofs are included in Appen-
dices. We note that this is arranged such that
our key ﬁndings 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 inﬁnite 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 inﬁnite 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 suﬀer from some kind of decoher-
ence, creating an inﬁnite 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 ﬁnite 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
ﬂuctuations of the probe observable, is deﬁned
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 ﬂuctu-
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 deﬁned 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-ﬂuctuations 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
inﬁnite 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 ﬁt 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, deﬁned analogously to the average
B
is the
bath Hilbert space dimension, and ﬁnally, 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 ﬂuctuation-
dissipation relation, which relates the time-
ﬂuctuations in the steady-state with the decay
rate after a quantum quench. The ratio between
ﬂuctuations 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
t
+ O
DE
, (12)
where O
DE
is the diagonal ensemble of the observ-
able O, which we ﬁnd 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 speciﬁc case of interest
here, Eq. (12) becomes,
hσ
z
(t)i = e
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 inﬁnite temperature
ﬂuctuation-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 inﬁnite-
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 eﬀective 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 eﬀects 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, deﬁned with an
energy dependent coupling strength, is the ther-
mal average decay rate over the energy width de-
ﬁned 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 eﬀect 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 ﬁnite 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 diﬀers 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 ﬂuctuation-dissipation theorem, relating
the time-ﬂuctuations 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 ﬂuctuations
δ
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
conﬁrms 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
)+
N1
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 ﬁts 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 conﬁrmation 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
inﬁnite temperature case, β = 0, and thus conﬁrms 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 deﬁned
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 ﬁnite 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-oﬀ
(ρ
B
)
αα
e
βE
α
to the bath state occupation.
For ﬁnite 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 inﬁnite 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 inﬁnite 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 inﬁnite temperature case, β = 0, and thus con-
ﬁrms the exponential scaling of χ(N)D(E)
1
N
B
. In
this case, we observe an exponential scaling for only
high temperatures, and conﬁrm 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 deﬁned
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
N1
, 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
N1
, 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-
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 diﬀers 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
N1
, 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 ﬁt),
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 [4446], can be
used to obtain arbitrary correlation functions
hc
µ
(α)c
ν
(α) · · · i
V
, where · · 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 brieﬂy 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 deﬁned as
hc
µ
(α)c
ν
(β)i
V
=Λ(µ, α)δ
αβ
δ
µν
(20a)
Λ(µ, α) :=
ω
0
Γ
(E
α
E
µ
)
2
2
(20b)
with Γ =
πg
2
N ω
0
. Λ
(2)
(µ, ν) is deﬁned similarly to
Eq. (20b), with Γ . 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 eﬀective 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 diﬀerent 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 eﬀective 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
diﬀerent 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 satisﬁed for such local observables.
The key assumption is related to the behaviour of
matrix elements, which must have a well deﬁned
average, that does not vary pathologically with
energy. A more formal deﬁnition 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-
ﬁne 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 fulﬁl 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 ﬁnd observ-
ables that do not fulﬁl 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 fulﬁlling
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-
ﬁned 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 deﬁned 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-
ﬁne
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
t
X
µν
µ6=ν
X
α
[O
αα
]
µ
w
αα
e
i(E
µ
E
ν
)t
Λ(µ, α)Λ(ν, α),
(34)
where to obtain the ﬁrst 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 deﬁning ˜µ = µ α, to obtain,
O(t) = hO(t)i
0
e
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 ﬁnd,
O(t) = hO(t)i
0
e
t
X
α
[