Modelling equilibration of local many-body quantum systems
by random graph ensembles
Daniel Nickelsen
1,2,3
and Michael Kastner
1,2
1
National Institute for Theoretical Physics (NITheP), Stellenbosch 7600, South Africa
2
Institute of Theoretical Physics, Department of Physics, University of Stellenbosch, Stellenbosch 7600, South Africa
3
African Institute for Mathematical Sciences, Muizenberg, Cape Town, South Africa
We introduce structured random matrix en-
sembles, constructed to model many-body
quantum systems with local interactions.
These ensembles are employed to study equili-
bration of isolated many-body quantum sys-
tems, showing that rather complex matrix
structures, well beyond Wigner’s full or
banded random matrices, are required to faith-
fully model equilibration times. Viewing the
random matrices as connectivities of graphs,
we analyse the resulting network of classical
oscillators in Hilbert space with tools from net-
work theory. One of these tools, called the
maximum flow value, is found to be an excel-
lent proxy for equilibration times. Since maxi-
mum flow values are less expensive to compute,
they give access to approximate equilibration
times for system sizes beyond those accessible
by exact diagonalisation.
1 Introduction
By means of typicality arguments, equilibration of iso-
lated quantum systems has been rigorously proved to
occur under rather general conditions [1–8]. It re-
mains however a key open problem to identify which
properties of quantum systems are responsible for en-
suring that equilibration timescales are physically re-
alistic [9–18]. As in virtually all other problems of
many-body quantum physics, the exponential growth
of Hilbert space dimension with system size renders
direct studies infeasible already for moderate system
sizes. To circumvent this difficulty in models of nuclei
of heavy atoms, Wigner proposed the use of random
matrices as Hamiltonians that statistically share cer-
tain properties of the nuclei [19–21]. One such choice,
the so-called Gaussian ensembles, is given by matri-
ces where all elements are normally distributed i.i.d.
numbers normalized to some measure [22]. Gaus-
sian random matrices are a particularly convenient
choice, as they allow for the analytic calculation of
ensemble-averaged properties like limiting distribu-
tions of eigenvalue statistics, i.e. Wigner’s semi-circle
Daniel Nickelsen: danielnickelsen@sun.ac.za
Michael Kastner: kastner@sun.ac.za
law [21]. Also, the level-spacing statistics of Gaussian
random matrices are known to follow the Wigner-
Dyson distribution [7, 22], which implies that small
spacings between eigenvalues of the Hamiltonian are
rare.
Gaussian ensembles of random matrices have also
been used to study equilibration and thermalisation of
isolated quantum systems [23–25]. Since the time evo-
lution of expectation values of observables is governed
by differences of energy eigenvalues, the level-spacing
statistics of random matrices influences the equilibra-
tion times [26]. The suppression of small spacings, as
in the Wigner-Dyson distribution, favours short equi-
libration times. Indeed, in Ref. [11] it has been shown
that equilibration under random-matrix Hamiltonians
typically happens extremely fast, much faster than
everyday experience or the analysis of models from
condensed matter theory suggests. This finding sug-
gests that Gaussian random matrices are unsuitable
for representing physically realistic quantum many-
body Hamiltonians.
Locality, in the sense of short-ranged interac-
tions between only few constituents in a spatially
extended many-body system, is believed to be a
key requirement for realistic equilibration dynamics.
Gaussian random matrices typically represent fully-
connected systems with distance-independent interac-
tion strengths between constituents. To incorporate
locality into random matrix ensembles, early works in
nuclear physics have introduced the two-body random
ensemble [27, 28]. This ensemble uses Gaussian ran-
dom matrices to model the two-body interactions in a
fermionic systems, and assembles these matrices into
the total Hamiltonian of a many-body system, which
turns out to be a sparse matrix. The two-body ran-
dom ensemble has also been used to study thermali-
sation of isolated fermionic systems [29–32]. General-
isations of the two-body random ensemble to random
k-body interactions have been studied for fermionic
as well as bosonic systems [33–35]. More generally,
random matrix ensembles that use full random ma-
trices for local interactions to assemble random total
Hamiltonians are referred to as embedded ensembles.
If one examines the structure of matrices from the
embedded ensembles, one finds that they are sparse
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 1
arXiv:1912.02043v4 [quant-ph] 25 May 2020
matrices and that they display a banded structure.
Similarly, the (non-random) Hamiltonian of a many-
body system with local interactions, when represented
in the tensor product basis with respect to which lo-
cality is defined, turns out to be a banded, sparse ma-
trix. In fact, already Wigner in his early attempts to
model nuclei of heavy atoms introduced banded ma-
trices that have nonzero elements only up to a certain
distance from the diagonal [19, 20]. In addition to
resembling Hamiltonians with local interactions, sev-
eral studies have shown that random band matrices
are also capable of modelling certain aspects impor-
tant for equilibration timescales, e.g. the transition
from Wigner-Dyson statistics for chaotic system to
Poisson statistics of regular systems, and the emer-
gence of many-body localisation [36–38]. Variants of
Wigner band matrices include power-law banded ran-
dom matrices [39, 40] where the standard deviation of
matrix elements decreases with the distance from the
main diagonal according to a power law, and sparse
random band matrices [41], which also account for the
possibility of having vanishing elements in the band.
However, while the mentioned ensembles of banded
random matrices may model certain aspects of physi-
cal systems more realistically, they come with the dis-
advantage that analytic calculations of their spectral
properties are far more intricate [42–44]. An overview
of applications of embedded and banded ensembles is
provided in Ref. [45]. In a different attempt to in-
clude aspects of locality into random Hamiltonians,
the authors of Ref. [46] disregarded the banded struc-
ture of local Hamiltonians, and instead investigated
only the effect of the sparseness of the matrix on the
equilibration dynamics.
In this paper we address in a systematic way the
interconnection between the locality of a many-body
Hamiltonian and the structural properties (banded-
ness, sparseness, etc.) of its matrix representation.
To this purpose, we interpret the Hamiltonian ma-
trix as the connectivity matrix of a weighted graph,
which allows us to apply concepts from graph and net-
work theory. As our first main result, we identify key
features of the graphs corresponding to local Hamil-
tonians, including degree distributions as well as a
node-dependent (“curved”) bandedness of the connec-
tivity matrix. As a second main result we identify
the matrix (or graph) properties that are essential to
faithfully reproduce the typical equilibration times of
a local many-body system. We find that a random
Hamiltonian with the correctly curved bandshape and
the correct, constant degree faithfully reproduces the
equilibration times of a local system. Abandoning ei-
ther the correct bandshape or constant degree, how-
ever, spoils the success. We also demonstrate that
Hamiltonians from such ensembles can be constructed
efficiently. Our third main result is that the maximum
flow, a tool from network theory, can serve as a proxy
for equilibration times. This means that, instead of
computing the equilibration time of a Hamiltonian, it
is sufficient—and numerically more efficient—to cal-
culate its maximum flow value, from which the equili-
bration time can be deduced. Combining the efficient
graph construction and the calculation of the maxi-
mum flow value, equilibration times of local Hamilto-
nians can be investigated more efficiently.
2 Equilibration
We choose the equilibration dynamics of a quan-
tum many-body system as our diagnostic tool with
which we analyse whether the locality of interactions
is successfully modelled by a random matrix ensemble.
Defining equilibration in an isolated, hence unitarily
evolving, quantum system requires some care, as we
discuss in the following. We consider the Schrödinger
time evolution
i∂
t
|ψ(t)i = H|ψ(t)i, (1)
generated by a time-independent Hamiltonian H. Un-
der the dynamics of (1), the infinite time average of
the expectation value of an observable O can be writ-
ten as
lim
T →∞
1
T
Z
T
0
hψ(t)|O|ψ(t)idt
=
X
j
|c
j
|
2
he
j
|O|e
j
i ≡ hOi
diag
, (2)
where |e
j
i are energy eigenstates and c
j
= he
j
|ψ
0
i are
the overlaps with the initial state. Unlike the micro-
canonical ensemble that uses equal weights, the so-
called diagonal ensemble average in the second line of
(2) uses the weights |c
j
|
2
. If the system equilibrates,
it is reasonable to expect that the expectation value
of the observable O equilibrates to the value hOi
diag
.
Through the dependence of the coefficients c
j
on ψ
0
,
the equilibrium value will in general depend on the
initial state.
In any system with a finite-dimensional Hilbert
space, unitary time evolution implies that hOi evolves
quasi-periodically in time. As a consequence, even at
arbitrarily late times, fluctuations around the equi-
librium value hOi
diag
occur. Because of these fluctua-
tions, a sensible definition of equilibration in isolated
macroscopic quantum systems requires a probabilis-
tic perspective. The largest part of the Hilbert space,
referred to as the equilibrium subspace, is spanned by
eigenstates |o
j
i of O that have eigenvalues o
j
close
to the equilibrium value, o
j
≈ hOi
diag
. The comple-
ment of this subspace is the nonequilibrium subspace,
spanned by eigenstates of O with eigenvalues that ap-
preciably differ from hOi
diag
. Equilibrium states are
then defined as states that have a large overlap with
the equilibrium subspace; see [11, 47] for details. As
a consequence of this definition, an equilibrium state
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 2
|ψi not only guarantees that hψ|O|ψi ≈ hOi
diag
, but
also that fluctuations around this expectation value
are small. For a many-body system with large Hilbert
space dimension, the dimension of the equilibrium
subspace is typically much larger than the dimension
of the nonequilibrium subspace. As a consequence, a
state drawn randomly from the Hilbert space accord-
ing to the Haar measure will most likely be an equilib-
rium state [11, 47]. Moreover, a nonequilibrium initial
state typically evolves towards equilibrium and stays
there for most times [3, 4]. Sporadic recurrences to
nonequilibrium occur, but with a probability that be-
comes vanishingly small with increasing system size.
Throughout this section, equilibration was dis-
cussed with reference to a specific observable O, and
also the (non)equilibrium subspaces have been de-
fined with respect to that observable. The reason for
this choice is related to the quasi-periodic behaviour
of a unitarily evolving quantum system on a finite-
dimensional Hilbert space. In such a system it is
always possible to find observables (i.e., Hermitian
operators) that either do not evolve at all (like the
projectors |e
j
ihe
j
| onto energy eigenspaces), or that
evolve periodically and hence do not equilibrate (like
|e
i
ihe
j
| + |e
j
ihe
i
|). The physical relevance of these
and similar observables is, however, questionable, as
they are nonlocal quantities (as defined in Sec. 4) and
usually inaccessible in experimental measurements of
many-body systems. Defining equilibration with re-
spect to specific observables is hence not a shortcom-
ing, but a desired feature.
3 Network picture of equilibration
The fact that equilibration in Sec. 2 has been defined
with respect to a specific observable suggests that, for
the physical problem at hand, the eigenbasis |o
j
i of
O may play a distinguished role. In this basis, the
Schrödinger equation (1) reads
˙x
j
(t) = −i
X
k
H
jk
x
k
(t), (3)
where
x
j
(t) = ho
j
|ψ(t)i, H
jk
= ho
j
|H|o
k
i. (4)
The eigenstates |o
j
i are assumed to be sorted in as-
cending order, o
j+1
≥ o
j
. By componentwise integra-
tion, Eq. (3) can be written as
e
iH
jj
t
x
j
(t) = x
j
(0) − i
X
k6=j
H
jk
Z
t
0
e
iH
jj
τ
x
k
(τ) dτ.
(5)
Taking the diagonal entries H
jj
as frequencies and the
off-diagonal entries H
jk
as coupling strengths between
classical oscillators, we can interpret this equation as
the description of a network of coupled oscillators with
amplitudes x
j
(t). These observations motivate to use
graph-theoretical tools for the analysis of the equili-
bration dynamics.
Using this idea of coupled oscillators, the network
picture of equilibration can be described as follows
[48]. Preparing the system in an initial state where
oscillators of equilibrium states have vanishing am-
plitudes x
j
(t), oscillations need to propagate through
the network in order to activate the equilibrium oscil-
lators. Once these oscillators are active, measurement
probabilities x
j
for equilibrium values o
j
≈ hOi
diag
are dominating and, according to the definition of
equilibrium in Sec. 2, the system has equilibrated.
The efficiency or speed of equilibration is expected to
depend on the network properties like degree distri-
bution, connectivity, etc. A more detailed motivation
of such a network picture of equilibration is given in
Appendix A. The main goal of the present paper is to
analyse the properties of such networks for physically
realistic quantum many-body systems and link those
network properties to the equilibration behaviour.
We consider a weighted graph consisting of N
nodes, where N is the dimension of the Hilbert space
of the quantum system. To each node j we assign the
classical complex degree of freedom x
j
defined in (4).
The weights of the graph edges are given by |H
jk
|.
The graph topology is described by the adjacency ma-
trix A with matrix elements
A
jk
=
(
0 if H
jk
= 0,
1 else.
(6)
We use
d
O
(j) :=
o
j
− hOi
diag
(7)
as a distance of the nodes to equilibrium, which will
serve as a metric on the graph: nodes j with d
O
(j) ≈ 0
are equilibrium nodes, and those with d
O
(j) 6≈ 0 are
nonequilibrium nodes.
In Ref. [48], this network picture of quantum dy-
namics has been used to derive a lower bound on the
timescale at which the system equilibrates. Remark-
ably, the bound was found to depend on the degrees
of locality of the Hamiltonian H as well as of the ob-
servable O. There exist various notions of locality,
and we will introduce two of them in Sec. 4, but we
expect all of them to have in common that the corre-
sponding network or graph is less strongly connected
when both H and O are local. Weaker and/or fewer
connections are then expected to lead to longer equili-
bration timescales, and this is indeed what was found
in Ref. [48].
There exists a zoo of measures of graph topology,
and our aim is to identify a measure that can serve
as a proxy for equilibration timescales. In view of
the above described intuitive picture of an oscilla-
tion of a nonequilibrium node propagating through
the network to an equilibrium node along all of the
existing pathways between the nodes, promising can-
didates are measures that take into account multiple
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 3
pathways. The node connectivity, defined as the min-
imum number of nodes that need to be removed in
order to split the graph into two disconnected sub-
graphs, takes multiple paths into consideration. A
closely related graph measure is the maximum flow
value f
max
, for which the weights |H
jk
| of the graph
are taken as capacities rendering the graph a flow net-
work, where f
max
is the total capacity between two
nodes. Visualizing the flow as a fluid entering the
network at node k and leaving it at node j, f
max
(k, j)
is the maximum amount of fluid per unit time that
can be routed through the network using those two
nodes.
Most algorithms use the max-flow min-cut theorem
to calculate f
max
. The theorem considers cuts of the
network into two disconnected parts. A capacity value
is assigned to each cut, which is the sum of all |H
uv
|
that have been removed for the cut. The minimum cut
is the one that is minimal in this capacity value. The
max-flow min-cut theorem then states that f
max
(k, j)
is identical to the capacity value of the minimum cut
separating nodes k and j. Intuitively, this means that
the maximum flow that can be routed through the
network is defined by the smallest bottleneck in ca-
pacity. A more formal definition of f
max
(k, j) can be
found in Ref. [49]. We will apply the maximum flow
measure to the network picture of quantum spin mod-
els in Sec. 5.3.
4 Spin model
There exist various notions of locality, and our goal
is to explore network features and exploit graph mea-
sures for quantum many-body systems with different
degrees of locality. To this aim, we introduce a rather
general two-parameter family of quantum spin chains,
where each of the two parameters d and n controls the
degree of locality according to one of those notions
of locality. The parameter d quantifies the degree
of locality according to the convention of, e.g., sta-
tistical physics or condensed matter physics, where
nearest-neighbour interactions correspond to d = 1,
next-nearest-neighbour interactions to d = 2, and so
on. The parameter n quantifies locality in the sense
of the number of lattice sites involved in an interac-
tion term, independently of the distance between the
sites. A noninteracting Hamiltonian corresponds to
n = 1, a Hamiltonian consisting of pair interaction
terms corresponds to n = 2, etc.
Consider a quantum spin chain on a one-dimension-
al lattice L consisting of L sites. (Not to be confused
with the nodes of the graph introduced in Sec. 3.) For
any subset l ⊆ L we denote by diam(l) the diameter
of l, i.e., the largest Euclidean distance between any
two sites in l. A general Hamiltonian for a spin-1/2
chain on L with locality degrees d and n is given by
H =
d
X
δ=1
X
χ(n)|δ
X
φ(n)
a(χ, φ)
n
Y
i=1
σ
φ
i
χ
i
, (8)
where χ(n)|δ denotes the set of all n-element subsets
l of L with the condition that diam(l) = δ. The
set φ(n) consists of all n-tuples of orientations x, y,
or z. σ
x
i
, σ
y
i
and σ
z
i
denote the x, y and z compo-
nents of Pauli spin operators acting on lattice site i.
The coupling constants a(χ, φ) are normal-distributed
i.i.d. numbers with zero mean and standard deviation
1/2. The random character of the couplings guar-
antees that accidental cancellations in the matrix el-
ements H
jk
occur with zero probability, and hence
the elements of the adjacency matrix (6) do not van-
ish by chance. Hamiltonians with parameter values
d = 1 and n = 2 correspond to pair-interactions where
each spin component of one spin interacts with all
spin components of only its nearest neighbours. For
d = L − 1 and n = L the Hamiltonian is the sum of
products σ
φ
1
1
σ
φ
2
2
···σ
φ
L
L
over all sets φ, in which case
the Hamiltonian is maximally nonlocal.
Our aim is to compare the equilibration timescales
of Hamiltonians of the type (8) for various values of
the parameters d and n. To make sure that such
a comparison focusses on the degree of locality, we
normalise all Hamiltonians such that kHk = 1 to
get rid of overall factors in H, as a Hamiltonian 2H
would equilibrate twice as fast as H for trivial rea-
sons. Unnormalised Hamiltonians usually scale ex-
tensively or superextensively with the lattice size L,
and our normalisation convention will therefore lead
to longer equilibration times compared to the unnor-
malised counterpart. We will comment on implica-
tions of the normalisation at the end of Sec. 5.2.
The magnetisation is a natural observable in a
quantum spin system. We study equilibration with
respect to a weighted variant of the magnetisation in
z-direction,
M =
1
L
L
X
i=1
a
m
(i) σ
z
i
, (9)
with real weights a
m
. The eigenstates |m
j
i of M
are products of eigenstates of σ
z
i
. For the case
of homogeneous weights, a
m
(i) ≡ 1 for all i, the
eigenvalues m
i
of M take on values from the set
{−L, L + 2, . . . , L − 2, L}/L. Most of these eigen-
values occur with high degeneracies, and hence the
eigenstates |m
i
i cannot be uniquely sorted according
to their eigenvalues. To avoid this, we introduce small
variations by choosing the a
m
as normal-distributed
random numbers with mean 1 and standard devia-
tion 1/10. A unique order of eigenstates brings out
the graph properties more clearly and demonstrates
that our analysis works for generic observables. Lift-
ing the degeneracy, however, is not a requirement for
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 4
our analysis, as is demonstrated in Appendix E. Due
to the independent, anisotropic and inhomogeneous
coupling, we expect the equilibrium value hMi
diag
to
be close to zero, which is confirmed by exact diago-
nalisation. Hence, equilibrium eigenstates are at the
centre and nonequilibrium eigenstates at the edges of
the spectrum.
5 Results
5.1 Graph properties
In this section we mostly study the properties of the
(unweighted) graph defined through the adjacency
matrix A (6) that is associated with the Hamiltonian
(8) when represented in the eigenbasis of the magneti-
sation M (9). We focus in particular on the depen-
dence of the graph properties on the locality parame-
ters n and d introduced in (8).
The degree distribution is a prominent graph prop-
erty. The degree %(j) =
P
k6=j
A
jk
of a node j counts
the number of edges of this node. The degree distri-
bution is the probability distribution of the degrees
over all nodes j. Scale-free and small-world networks,
for instance, have power-law distributed degrees, and
the degrees of Bernoulli random graphs follow a Pois-
son distribution. Here we find that, for any given n
and d, all nodes of the graph have identical degrees.
Such a graph is called a regular graph. By examining
constant degrees % for system sizes up to L = 10 and
all possible combinations of locality parameters n and
d (and up to L = 16 for small values of n and d), we
conjecture that
%(L, n, d) =
n−1
X
q=0
L − q
d + 1
q + 1
d
q
. (10)
See Appendix B for details.
For a first impression of the effect of locality on the
adjacency matrix A, Fig. 1 shows graphical represen-
tations of the matrix elements A
jk
for different choices
of the parameters n and d. The nearest-neighbour
pair-interaction Hamiltonian [n = 2 and d = 1 in (8),
top left in Fig. 1] is a sparse matrix with a rather low
density of nonzero elements. The maximally nonlocal
Hamiltonian (n = L and d = L − 1, bottom right in
Fig. 1), on the other hand, is a dense matrix, which
may be modelled by full random matrices of the Gaus-
sian orthogonal ensemble. Other choices of n and d
create matrices of intermediate sparsity.
The plots in Fig. 1 reveal that, in addition to be-
ing sparse, the adjacency matrices of sufficiently lo-
cal systems are banded, in the sense that, outside a
certain region around the main diagonal, all matrix
elements are zero. Such a banded structure agrees
with Wigner’s intuition that banded random matri-
ces may be more suitable to model physical systems
[19, 20]. Also, a banded structure is consistent with
Figure 1: Plots of the adjacency matrices A corresponding
to the Hamiltonian (8) in the eigenbasis of the magnetisation
(9) for a chain of L = 8 spins and all possible combinations of
values for n and d. Black indicates nonzero matrix elements,
white indicates zero elements.
a recent theorem by Arad, Kuwahara, and Landau
[50], which affirms that a local observable represented
in the eigenbasis of a local Hamiltonian is a banded
matrix, possibly decorated with exponentially sup-
pressed nonzero elements outside the band. The rea-
soning that led to this result can be reversed to show
that a local Hamiltonian is banded in the eigenbasis
of the observable, again with the possibility of ex-
ponentially small matrix elements outside the bands
[48]. However, for the choice of Hamiltonian and ob-
servable considered in the present paper, we find that
the matrix elements H
jk
are strictly zero outside the
band.
By inspection of the upper rows of Fig. 1, we ob-
serve that the bands of nonzero matrix elements are
curved, in the sense that in the centre of the matrix
nonzero elements occur further away from the diago-
nal compared to the edges of the matrix. We call this
a node-dependent bandwidth. Since equilibrium states
are located towards the centre of the matrix, this im-
plies that a nonequilibrium node j is connected only
to nodes with a similar index i ≈ j, whereas equi-
librium nodes may also be connected to nodes whose
index i differs more substantially from j. For the ad-
jacency matrix of the spin Hamiltonian (8), the node-
dependent bandwidth can be calculated analytically.
We first derive such a result for the magnetisation (9)
with homogeneous weights a
m
(i) ≡ 1, and then gen-
eralise the argument to arbitrary observables.
Writing the Hamiltonian in terms of ladder oper-
ators σ
±
i
=
1
2
(σ
x
i
± iσ
y
i
), n-spin terms translate into
products of n ladder operators, which allow for at
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 5
Figure 2: Bandwidth functions b(j) and b
m
(j) compared to the actual distance
˜
b(j) and
˜
b
m
(j) of the nonzero entry furthest
from the diagonal of adjacency matrices A
jk
from (6) derived from Hamiltonian in (8) for n = 2 and d = 1 (left) and n = 3
and d = 2 (right). The bandwidth b(j) is defined in (12) and applies to the case a
m
(i) ≡ 1 for the magnetisation M in (9)
as observable. The bandwidth b
m
(j) follows from the general formula (14), here applied to the randomised M using normal-
distributed a
m
(i) as defined in (9). The bandwidths
˜
b(j) and
˜
b
m
(j) follow from the definition in (13) for non-randomised and
randomised M respectively.
most n spin flips and can change the magnetisation
by at most ∆m = 2n/L. We define the bandwidth
b(j) = max
k:o
k
≤o
j
−∆m
|k − j| (11)
of node j as the largest distance to any of the nodes
k for which o
k
≤ o
j
−∆m. For the magnetisation (9)
with a
m
(i) ≡ 1 there exist
L
q
states with q up-spins,
hence we can write
b(j) =
n
X
q=0
L
j + q
, (12)
where we assumed for simplicity that node j is the
first node of a block of constant magnetisation. Note
that b depends on the locality parameter n, but not
on d. Since H
jk
is typically not densely filled inside
the band, the actual distance of the nonzero entry
furthest from the diagonal,
˜
b(j) = max
k:A
j−k,j+k
6=0
|k − j|, (13)
may be smaller than the bandwidth,
˜
b ≤ b.
Noting that
L
q
is the density of eigenstates for
magnetisation m = (2q−L)/L, we can generalize (12)
to general observables O,
b
o
(j) =
Z
o
j
+∆o
o
j
g(o) do = G(o
j
+ ∆o) − G(o
j
). (14)
Here, g(o) is the density of eigenstates O, and G(o)
the integrated density of states, defined as the number
of eigenstates |o
j
i of O for which o
j
≤ o. By
∆o = max
ψ
hψ|H
†
OH − O|ψi (15)
we denote the Hamiltonian’s circle of influence in the
spectrum of the observable, which can be determined
by maximising over the eigenstates |o
j
i of O. Equa-
tion (14), being general, also applies to the magneti-
sation (9) with normally distributed weights a
m
(i)
that we are interested in. Figure 2 illustrates how
the bandwidths b(j) and b
m
(j) compare with the ac-
tual structure of the adjacency matrix A
jk
of the spin
model (8) when using the magnetisation (9) as ob-
servable. Overall, we find that the locality parameter
n determines the bandwidth of the adjacency matrix,
whereas n and d determine the density of nonzero el-
ements inside the band.
5.2 Comparison of random graph ensembles
In the previous section we have studied certain graph
features of adjacency matrices, in particular their de-
gree distributions % and node-dependent bandwidths
b, and investigated the dependence of these features
on the locality parameters n and d of the Hamilto-
nian. Now we reverse the logic by constructing ran-
dom graphs with a certain b and %, and investigate to
what extend the corresponding random matrix ensem-
bles reproduce the equilibration dynamics of Hamilto-
nians with locality parameters n and d. In the follow-
ing we sketch the main features of the random graph
ensembles considered in this paper. Full descriptions
on how the ensembles are constructed are given in
Appendices C and D.1.
EXH: Exact Hamiltonian H
jk
of the Heisenberg
model (8) with random couplings drawn from a
normal distribution with zero mean and standard
deviation 1/2.
EXA: Based on the exact adjacency matrix A of
the Hamiltonian (8), random Hamiltonians are
created by replacing the nonzero matrix elements
of A with normally distributed i.i.d. numbers.
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 6
EXH EXA BRF
BVF BRC REG
Figure 3: Plots of the adjacency matrices for the various
random graph ensembles for L = 8 spins. Note that EXH
and EXA are identical since the resulting Hamiltonians dif-
fer only in the statistics of their matrix elements, not their
adjacencies.
BRF: Randomly generated banded regular graphs
with constant degree %(L, n, d) and fixed
functional bandwidths b
m
(j). The functional
bandwidth (14) of the exact Hamiltonian (8) with
a given n and d is used.
BVF: Randomly generated banded graphs with
variable degree and fixed functional bandwidth
b
m
(j) using (14). The average degree is set to
%(L, n, d), and the functional bandwidth matches
that of the Hamiltonian (8) with a given n and
d.
BRC: Banded regular graphs with constant (node-
independent) bandwidth b ≡ max
j
b
m
(j) using
(14) are generated randomly but not uniformly.
The maximum bandwidth and the constant de-
gree match those of the original Hamiltonian (8).
REG: Regular graphs with constant degree
%(L, n, d), sampled uniformly. This ensemble is
similar to the approach in [46].
In Fig. 3 we show adjacency matrices for these six
random graph ensembles. The graphs of all the en-
sembles have the same total number of edges, i.e.
all resulting Hamiltonians have the same number of
nonzero matrix elements. Apart from the EXH en-
semble, the weights of all graph edges are i.i.d. num-
bers sampled from normal distributions matching the
statistics of the exact H
jk
. Details of the weight dis-
tributions are given in Appendix D.2.
We interpret each realisation of a random graph
drawn from one of the above ensembles as the matrix
representation H
jk
of a Hamiltonian in the eigenbasis
of the randomised magnetisation M in (9). Our aim
is to study the equilibration timescales of the Hamil-
tonians and assess which of the ensembles lead to
timescales that faithfully reproduce those of Hamil-
tonians with the corresponding degrees of localities n
and d. As discussed in Sec. 4, timescales are made
more comparable by normalising each H
jk
such that
kHk = 1. Since equilibrium values of M are expected
to be close to zero, we choose the initial state as far
away from equilibrium as possible, i.e. we select the
node k = N that corresponds to a maximum mag-
netisation.
Using the function expm_multiply from the python
package scipy.sparse.linalg, which is designed for
efficient propagation of initial vectors under the linear
time evolution of sparse matrices, we compute the ex-
act time evolution for samples of Hamiltonians drawn
from each of the above graph ensembles. For each
time evolution we measure the equilibration time T
eq
,
defined as the time at which both hOi and
O
2
for
the first time reach their diagonal averages with an
allowed margin of error of 10% of their respective ini-
tial expectation values [48]. This definition of T
eq
is
in the spirit of the definition of equilibration in terms
of equilibrium subspaces in Sec. 2.
Figure 4 shows T
eq
for various system sizes. Based
on these plots, we make the following observations:
(a) The exact ensemble EXH shows an approxi-
mately linear increase of T
eq
with the system size
L for all three sets of parameters shown in Fig. 4.
(b) The equilibration times obtained from the REG
and BRC ensembles are substantially smaller
than those of EXH and do not (or hardly) in-
crease with L. This disqualifies REG and BRC
as unsuitable for reproducing the equilibration
behaviour of local Hamiltonians.
(c) The BVF ensemble overestimates the increase of
T
eq
with L, showing slower equilibration than
EXH, at least for the larger system sizes. We
speculate that variations in the degrees of nodes
cause a bottleneck effect, where poorly connected
nodes do not permit efficient propagation of os-
cillations towards equilibrium nodes.
(d) Unsurprisingly, the equilibration times of EXH
agree best with those of EXA, the ensemble
which has the exact same adjacency matrix.
Their Hamiltonian matrix elements, however, dif-
fer, even though they are statistically equivalent
if taken to be i.i.d. numbers. The differences in
T
eq
can be attributed to correlations between the
matrix elements H
jk
that have not correctly been
accounted for in EXA. This observation is in line
with recent studies on the role of correlations be-
tween the random matrix elements of the Hamil-
tonian and the observable [51–53], a topic that
has been discussed also as an extension to the
eigenstate thermalisation hypothesis [54–56]. To
fully focus on the effects of graph topology, it
may make sense to compare equilibration times
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 7
Figure 4: Dependence of the equilibration times T
eq
on the system size L for the six random graph ensembles introduced
in Sec. 5.2. Averages and standard errors of T
eq
are computed using sample sizes decreasing exponentially as 2
18−L
. The
smaller sample sizes are justified due to the self-averaging properties of the larger matrices. Different disorder realisations of
the randomised magnetisation (9) are used for all random graphs. Our computational resources restricted us to system sizes
of only L = 10 for the computationally more demanding EXH ensemble. The range of the ordinate focusses on relevant values
for T
eq
for clearer plots.
of the various ensembles to EXA instead of EXH,
so that different topologies but identical statistics
of matrix elements H
jk
are compared.
(e) Among the ensembles with random adjacency
matrices, the BRF ensemble is the only ensemble
reproducing the linear scaling of T
eq
with L ob-
served for the ensembles EXH and EXA, and as
such best suited to model the equilibration dy-
namics for the family of Hamiltonians considered
in (8). However, on the basis of the available
data, one may speculate that the quantitative
agreement deteriorates with decreasing locality.
We conclude that degree distribution and banded-
ness are essential features of the adjacency matrices
of local Hamiltonians that need to be taken into ac-
count in a random matrix approach. The shape of the
node-dependent bandwidth also has a role to play, but
appears to be somewhat less crucial. Since it is com-
putationally costly to calculate the exact adjacency
matrix required for the EXA ensemble, resorting to
BRF as the next-best, and computationally less ex-
pensive, ensemble is a potential avenue for studying
equilibration times of larger systems in a random ma-
trix setting. We do however not claim that the BRF
ensemble with its constant degree is always the best
choice to model equilibration of local systems. Many
other natural Hamiltonians may not be of constant
degree (e.g. particular choices of the a(χ, φ) in (8)) –
studying the degree distribution of these systems and
how it effects their equilibration timescales should be
interesting for future research.
A comment is in order on the increase of T
eq
with
system size L that is observed in EXH as well as in
several of the approximate random graph ensembles.
On physical grounds, one expects a (statistically)
homogeneous local system with homogeneous initial
conditions to equilibrate on an approximately con-
stant (system-size independent) timescale: because of
the homogeneity, no currents across macroscopic dis-
tances are required to reach equilibrium, and local
equilibration is sufficient to reach global equilibrium.
This reasoning applies in the conventional setting of
an extensive local Hamiltonian. Our convention of
normalising H such that kHk = 1, however, intro-
duces an extra factor of 1/kHk in the time-evolution
operator, resulting in the observed approximately lin-
ear increase of T
eq
with L.
5.3 Maximum flow vs. equilibration time
To determine the equilibration time for a given Hamil-
tonian, one needs to integrate the time-dependent
Schrödinger equation up to sufficiently late times.
This is computationally costly and is feasible only for
small system sizes. To circumvent this difficulty, we
propose to calculate instead the maximum flow value
f
max
, defined in Sec. 3, of the graph defined by H
jk
,
which is numerically less expensive. We show that
T
eq
and f
max
are strongly correlated, and by calculat-
ing the latter one can estimate the equilibration time
to fairly good accuracy. Extrapolating this finding to
larger system sizes, estimates of equilibration times of
larger systems can be obtained.
We use the preflow_push algorithm of the python
package NetworkX to calculate f
max
(k, j) for various
realisations of the spin Hamiltonian (8). In Fig. 5
(left) we show a scatter plot of the pairs (f
max
, T
eq
)
for various values of n and d, and in Fig. 5 (right) the
corresponding averages and standard deviations. The
plots show strong anticorrelations between f
max
and
T
eq
: shorter equilibration times are observed when
the maximum flow value is larger, and vice versa.
Small imperfections in the anticorrelations may point
towards limitations of the usage of f
max
as a proxy for
T
eq
, or they may simply be caused by shortcomings
of our numerical method for computing equilibration
times.
f
max
and T
eq
scale differently with system size L.
To use f
max
as an estimate for T
eq
for larger sys-
tems, we fix n and d and calculate f
max
and T
eq
for
random samples of Hamiltonians for system sizes be-
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 8
Figure 5: Left: Scatter plot of the maximum flow value f
max
versus equilibration time T
eq
for L = 10 spins, calculated for all
possible combinations of n = 2, . . . , 10 and d = n − 1, . . . , 9. Different symbols correspond to different values for n, different
colours to different values of d (see both legends). For each combination of n and d, 2
10−n
+ 1 instances of H
jk
from (8) are
used to obtain data sets of pairs (T
eq
, f
max
). Right: Averages and standard errors of the data in the left plot, separately for
each combination of n and d (standard errors are smaller than the symbol sizes).
Figure 6: Sample averages of maximum flow values f
max
versus equilibration times T
eq
for random Hamiltonians from the BRF
ensemble. Averages and standard errors of 2
14−L
samples are shown, where the standard errors for f
max
are smaller than the
symbol sizes. The straight lines are linear fits to the data.
tween 4 and 10. In Fig. 6 we show results for random
Hamiltonians from the BRF ensembles corresponding
to (n = 3, d = 2) and (n = 2, d = L − 1). The data
suggest a linear relation between T
eq
and f
max
. A lin-
ear fit to the data, shown in the plot, is used in Sec. 5.4
to estimate equilibration times of larger systems.
5.4 Equilibration times of larger systems
Based on the results of Secs. 5.2 and 5.3, we are now
in the position to sample Hamiltonians from the BRF
ensemble for different locality parameters, compute
the maximum flow values of the corresponding graphs,
and use the linear extrapolations of Fig. 6 to deduce
T
eq
for system sizes up to L = 17. While these are
still small systems, they are larger than those that can
be treated directly with the computational resources
at our disposal. In Fig. 7 we plot the thus obtained
results for system sizes between 4 and 17, and for
various combinations of locality parameters n and d.
As expected, T
eq
increases with L, and larger values
of n entail faster equilibration.
The limitation to L 6 17 is due to the calculation
of f
max
, which has a complexity of O(
√
%N
5/2
) when
using the highest-label preflow-push algorithm imple-
mented in NetworkX. The algorithm used to gener-
ate the random graphs in the BRF ensemble (see Ap-
pendix D.1) has an estimated complexity of O(%N),
with potential for improvement.
6 Conclusions
We introduced and studied various random graph en-
sembles, with the aim of modelling equilibration in a
spin system with a given degree of locality. Interpret-
ing the Hamiltonian matrix elements as the weights of
a graph in Hilbert space, we employed concepts and
tools from classical network theory to study quantum
mechanical time evolution. The rationale of our anal-
ysis proceeds in three steps:
(a) Working with a local Hamilton H and a local
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 9
Figure 7: Equilibration times T
eq
obtained from maximum
flow values f
max
of random graphs of the BRF ensemble for
various system sizes L. Different disorder realisations of the
randomised magnetisation (9) are used for all random graphs.
The relation between T
eq
and f
max
is given by linear fits as
shown in Fig. 6. Displayed are averages and standard errors
of max(4, 2
18−L
) samples for each L.
observable O implies the existence of a basis in
which H and O are simultaneously sparse. A nu-
merical study of quantum spin models with ran-
dom couplings revealed further manifestations of
locality in the matrix structure, in particular a
locality-dependent constant degree (10) and a
node-dependent bandwidth (12).
(b) Based on these observations we introduce the
random graph ensembles EXA, BRF, BRC, BVF,
and REG, which incorporate degree distribution
and bandwidth to a certain extent. We numer-
ically study samples of Hamiltonians from the
various ensembles, finding that the BRF ensem-
ble faithfully captures the dependence of equili-
bration times on locality, while having compu-
tational advantages over the ensemble EXH of
exact random Hamiltonians.
(c) To avoid the costly numerical computation of
equilibration times, we showed that the max-
imum flow value f
max
between nonequilibrium
and equilibrium nodes of a graph is strongly cor-
related with the equilibration time T
eq
of the cor-
responding Hamiltonian. Calculating f
max
and
inferring T
eq
is computationally less costly and
gives access to equilibration times for system sizes
that cannot be reached by direct numerical inte-
gration.
Even though we started out from the assumption
that the random graphs represent the Hamiltonian in
the eigenbasis of the observable, the construction of
the random graph ensemble is not specific to a sin-
gle observable: only the observable’s density of states
g(o) and the Hamiltonian’s circle of influence ∆o (15)
enter the construction, and all observables that share
these features can be modelled by the same random
graph ensemble. The effect of g(o) and ∆o on equili-
bration times for large system sizes can then be stud-
ied using f
max
of the generated graphs, without deriv-
ing exact Hamiltonians, spectral properties, or solving
the Schrödinger equation.
Presently, the modelling of local systems by random
graph ensembles, and also the estimation of equili-
bration times through maximum flow values, led to
only a moderate increase of the system sizes we can
deal with, but there is plenty of potential for further
improvement. The method for constructing banded
random adjacency matrices with constant degree de-
scribed in Appendix D.1 is certainly far from optimal,
and also does not sample uniformly, which may affect
the accuracy of the modelling. Also, graph measures
other than the maximum flow value may be more
suitable and are worth being investigated. We have
looked into a few others, including the node connec-
tivity and the shortest path length, but these turned
out to be less strongly correlated with T
eq
.
Finally, we believe that the random graph en-
sembles introduced in this paper, as well as the
network-theoretic tools employed to analyse the ran-
dom graphs, have potential for the analysis of local
quantum systems well beyond the specific question of
equilibration addressed here. It should be interesting
to study the effects of locality on the level spacing
statistics, a topic that we are planning to report on
in future work.
Acknowledgements
The authors gratefully acknowledge helpful comments
on the manuscript from Fausto Borgonovi. M. K. ac-
knowledges financial support from the National Re-
search Foundation of South Africa through the In-
centive Funding and the Competitive Programme for
Rated Researchers.
A More on the network perspective
In Ref. [48] it has been argued that the quantity
Λ
jk
(t) =
∂x
j
(t)
∂x
k
(0)
=
e
−iHt
jk
=
∞
X
q=0
− ikHkt
q
q!
(H
q
)
jk
kHk
q
,
(16)
can be used to measure the speed of equilibration.
Λ
jk
(t) quantifies the influence a nonequilibrium node
k at time zero has on an equilibrium node j at time
t. Typically, Λ
jk
(t) is small for small values of t, and
increases with time until it saturates. A large value
of Λ
jk
indicates that a change in x
k
affects x
j
after a
time t. If the graph distance between nodes k and j
of the network is maximal, then any initial state can
excite the equilibrium node x
j
, as all other nodes are
closer to j. Therefore, when Λ
jk
(t) saturates close to
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 10
its maximal values, we can say the system has equili-
brated.
It is insightful to have a closer look at the expansion
(16) in powers of H. The matrix components of these
powers are of the form
(H
q
)
jk
=
X
i
1
,...,i
q−1
H
ji
1
H
i
1
i
2
···H
i
q−2
i
q−1
H
i
q−1
i
k
=
X
ω
jk
(q)
Y
κ∈ω
jk
(q)
H
κ
=
X
ω
jk
(q)
C (ω
jk
(q)), (17)
where the sums in the second line run over all walks
ω
jk
(q) that connect nodes j and k using q steps, e.g.
ω
25
(4) ∈ {{(24), (43), (34), (45)}, . . .}. Walks are al-
lowed to visit the same nodes multiple times. To each
walk we assign the corresponding coupling chain
C (ω
jk
(q)) =
Y
κ∈ω
jk
(q)
H
κ
, (18)
and hence Eq. (17) expresses the qth power of H as
the sum over all possible coupling chains of length
q. For sparse matrices H, all small-q coupling chains
between far-apart nodes may be zero. We can then
interpret the expansion (16) as the weighted sum of
all coupling chains C (ω
jk
(q)) between nodes j and k,
with q-dependent weights w
q
= (−ikHkt)
q
/q!. For
t = 0, nodes j and k are uncoupled. For small times
t < 1/kHk, |w
q
| peaks for couplings of short length
(q ≈ 2), and (H
q
)
jk
contributes to the sum only for
small values of q. Once t > 1/kHk, peak values of |w
q
|
shift to longer coupling chains (q 1), and (H
q
)
jk
with larger q contribute to the sum. Since the statis-
tics of the (H
q
)
jk
saturates for sufficiently large q,
and due to the alternating signs of (−i)
q
, the value
of Λ
jk
(t) oscillates around a stable average for suf-
ficiently large t. This is the saturation mentioned
above, which indicates that the system has equili-
brated.
The inverse of the norm of the Hamiltonian, 1/kHk,
therefore gives an indication on the equilibration
timescale, but it does not capture the influence of lo-
cality. Locality enters through the vanishing coupling
chains of small lengths q. For the dense Hamilto-
nian of a nonlocal system, the weights w
q
on short
coupling lengths q ≈ 2 for small t are already suffi-
cient for the fluctuations of Λ
jk
(t) to saturate. For
the sparse Hamiltonian of a local system, however,
coupling lengths are longer and increase with system
size, and t needs to be sufficiently large in order for
the w
q
to give a nonnegligible weight to these coupling
lengths, resulting in longer equilibration times.
In Fig. 8 we illustrate how the two factors w
q
and
(H
q
)
jk
that enter the expansion (16) depend on q
for the Hamiltonian (8) with n-spin interactions over
a maximum range of d lattice sites. For nearest-
neighbour pair-interactions (n = 2, d = 1) and L = 8
spins, coupling chains of length q ≥ 2 contribute. For
a larger system of L = 12 spins only coupling chains
of length q ≥ 4 contribute, for which w
q
is negligi-
ble for t ≈ kHk. For less local choices of parameters
(n = 5, d = 6), |w
q
| already peaks at q = 2, bringing
the system close to equilibrium for t ≈ kHk.
B Constant degree
As a consequence of statistical isotropy and homo-
geneity of the Hamiltonian, all nodes of a given adja-
cency graph have a constant degree %. By numerically
determining % for system sizes up to L = 16 and var-
ious combinations of locality parameters n and d, we
conjectured that
%(L, n, d) = L
n−1
X
q=0
d
q
− C(n, d), (19)
where the first term depends linearly on L, and the
not yet determined second term C(n, d) is indepen-
dent of L. By studying C for various fixed values of
n we were able to further conjecture that
C(2, d) =
d−1
X
i=0
i
X
j=1
1
, (20a)
C(3, d) =
d−1
X
i=0
i
X
j=0
1 +
j
X
k=1
2
, (20b)
C(4, d) =
d−1
X
i=0
i
X
j=0
1 +
j
X
k=1
2 +
k
X
l=2
3
, (20c)
and so on. Extrapolating the observed pattern to ar-
bitrary n we obtain
C(n, d) =
d−1
X
j
1
=0
j
1
X
j
2
=0
1 + 2
d−1
X
j
1
=0
j
1
X
j
2
=0
j
2
X
j
3
=1
1 + 3
d−1
X
j
1
=0
j
1
X
j
2
=0
j
2
X
j
3
=1
j
3
X
j
4
=2
1 + ··· + (n −1)
d−1
X
j
1
=0
j
1
X
j
2
=0
j
2
X
j
3
=1
···
j
n−1
X
j
n
=n−2
1
=
d−1
X
j
1
=0
j
1
X
j
2
=0
1 +
d−1
X
j
1
=1
j
1
X
j
2
=1
j
2
X
j
3
=1
2 +
d−1
X
j
1
=2
j
1
X
j
2
=2
j
2
X
j
3
=2
j
3
X
j
4
=2
3 + ··· +
d−1
X
j
1
=n−2
j
1
X
j
2
=n−2
j
2
X
j
3
=n−2
···
j
n−1
X
j
n
=n−2
(n − 1)
=
X
0≤j
2
≤j
1
≤d−1
1 +
X
1≤j
3
≤j
2
≤j
1
≤d−1
2 +
X
2≤j
4
≤j
3
≤j
2
≤j
1
≤d−1
3 + ··· +
X
n−2≤j
n
≤j
n−1
≤···≤j
2
≤j
1
≤d−1
(n − 1)
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 11
Figure 8: Illustration of the terms in (16) with w
q
= (−it)
q
/q!, and H from the EXA ensemble introduced in Sec. 5.2 and
normalised such that kHk = 1. As indicated in the plots, the data are for different numbers of spins L, times t, and locality
parameters n and d as defined in Sec. 4. The (H
q
)
jk
and the |w
q
| are divided by their respective maxima for clarity.
=
n−1
X
q=1
q
X
q−1≤j
q+1
≤j
q
≤···≤j
1
≤d−1
1
=
n−1
X
q=1
q
d + 1
d − q
=
n−1
X
q=1
q
d + 1
q + 1
. (21)
In the last step of this calculation we used that
a+b−1
b−1
is the number of distinct ways to draw a
times from b options with replacement where the or-
der does not matter, such that with a = q + 1 and
b = d − 1 − (q − 1) + 1 we get the number of terms
in each sum. Inserting this expression into (19) we
obtain
%(L, n, d) =
n−1
X
q=0
L
d
q
− q
d + 1
q + 1
, (22)
which agrees with Eq. (10). We tested (10) for all
combinations of locality parameters n and d for sys-
tems of sizes up to L = 10, and also for system sizes
up to L = 16 for smaller values of n and d.
C Construction of Hamiltonians
The Hamiltonian (8) can be written as
H =
X
χ(n) | (χ
n
−χ
1
≤d)
h
χ
1
χ
2
...χ
n
. (23)
with
h
χ
1
χ
2
...χ
n
=
(x,y,z)
X
φ
1
,φ
2
,...,φ
n
a
φ
1
...φ
n
χ
1
...χ
n
2
χ
i
−1
⊗ σ
φ
1
χ
1
⊗
2
χ
2
−χ
1
−1
⊗ σ
φ
2
χ
2
⊗ ··· ⊗
2
χ
n
−χ
n−1
−1
⊗ σ
φ
n
χ
n
⊗
2
L−χ
n
−1
. (24)
In python, we can make use of the package
scipy.sparse. We first create sparse matrices for
σ
x
, σ
y
, and σ
z
in the σ
z
-eigenbasis, using the coordi-
nate format coo. The Kronecker products ⊗ in (24)
are available as scipy.sparse.kron, for efficiency
reasons the option formate=’coo’ should again be
used. Using Knuth’s “algorithm L”, all permutations
of a sequence with n 1’s and L − n 0’s are efficiently
created and each permutation is used to construct
h
χ
1
χ
2
...χ
n
(24). To create all tuples φ(n), the func-
tion product from itertools is used. Summing over
all thus generated terms and inserting the random
numbers a
φ
1
...φ
n
χ
1
...χ
n
yields the full Hamiltonian from the
EXH ensemble defined in Sec. 5.2. On the personal
computer we used, this procedure works for system
sizes up to L = 13.
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 12
Sheet1
Page 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
e z
1
o - - o o - o o o o - -
7 0
2
x - o o o o o - - - - - o
7 0
3
- - o o - o o - o - o - - o
7 0
4
- x x o o o - o o - - - - - -
7 0
5
x x x x - - - - o - o - o - -
7 0
6
x x - x - - - o o - - - o o -
7 0
7
- x x x - - o o o - - o - - -
7 0
8
x x x - - - x - o - - o - o -
7 0
9
x - - x - x x -
4 7
10
x - x x x x x x
7 7
11
x - - - - - - -
1 7
12
- - x - x - - -
2 7
13
- - - - - - x x
2 7
14
x - - x x - -
3 7
15
x - - x - x
3 7
16
- - - - -
0 7
Figure 9: An illustration of a conflict in the progressive con-
struction of regular graphs with banded adjacency matrices
for L = 4. The open circles in the upper triangle show
positions of randomly placed entries inside the bandwidth,
entries outside the band are shaded in grey. Crosses in the
lower triangle mark positions of entries that are required in
order to create a symmetric matrix. The two columns on
the right show the number of entries (e) and the number of
empty places (z) for each row. Here, the constant degree
is % = 7, so each row will have to have 7 nonzero entries.
The described routine goes through all rows from top to bot-
tom, randomly selects entries in the upper triangle and adds
entries in the lower triangle to ensure a symmetric matrix.
After completing the 8th row, the last row has no entries yet
and 7 empty places left and should therefore be completed
first before proceeding with row 9. However, the 10th row
is complete already, and hence either symmetry or constant
degree has to be sacrificed. To resolve this conflict, one of
the rows 1, 4, 5, 6, 7, 8 will be shuffled and other entries be
updated accordingly until the 10 row has an empty place such
that the last row can be completed.
By specifying the data type of the sparse matrices
as boolean, dtype=’bool’, and omitting the random
variable a
φ
1
...φ
n
χ
1
...χ
n
, mask matrices can be generated effi-
ciently for up to L = 18 spins on a standard computer.
The mask matrices are used as adjacency matrices for
the EXA ensemble defined in Sec. 5.2.
D Random graph ensembles
D.1 Construction of adjacency matrices
For the random graph ensembles BRF, BVF, BRC,
and REG introduced in Sec. 5.2, we need to construct
random adjacency matrices with a given constant de-
gree %. For the full (non-banded) regular graphs
in REG, the function random_regular_graph in the
python package NetworkX generates such matrices for
a given degree % and number of nodes N uniformly
from all possible matrices. We are not aware of simi-
lar ready-to-use packages for banded random matrices
of a prescribed bandshape. Here we describe a simple,
non-optimised algorithm to generate banded random
matrices with constant degree. Note that, as opposed
to random_regular_graph, our algorithm does not
ensure a uniform sampling of matrices.
The challenge is that the adjacency matrices are
required to be symmetric. For that reason, going
through all rows from top to bottom and randomly
choosing a fixed number of nonzero elements within
a prescribed bandwidth cannot work, as the corre-
sponding columns must be filled concurrently, which
in turn means that rows further down are partly filled
as we go. One therefore must keep track of the num-
ber e of nonzero elements already placed in each un-
completed row, and choose only as many nonzero en-
tries as required to obtain the constant degree for that
row. Moreover, the zero entries of each completed row
must remain zero in the corresponding columns, re-
ducing the number z of entries in the uncompleted
rows that can be nonzero. Before completing each
row with the missing number of nonzero elements, all
later rows must be checked for available space. Once a
row has just enough space left as needed to obtain the
desired constant degree for that row, it must be com-
pleted and the number of nonzero entries and space
for the uncompleted rows must be updated. When
the row has been completed and no other rows need
to be completed first, the routine resumes with the
topmost uncompleted row.
The algorithm described so far will most likely run
into a conflict: a row with just enough space to ob-
Figure 10: Exponential functions exp(aL
b
) with fit param-
eters a and b, fitted to disorder averages of ∆o (12) (top)
and kHk (bottom) for 100 disorder samples. Both plots are
for the Hamiltonian (8) with parameters n = 1, d = 1, and
system size L = 8.
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 13
Figure 11: Histograms of diagonal (real-valued) and off-diagonal (real and imaginary parts) matrix elements H
jk
of the
Hamiltonian (8) with n = 2 and d = 1 on L = 8 lattice sites. Statistics of matrix elements is obtained by generating 100
samples of Hamiltonians. Mean and standard deviations are estimated and the resulting normal distributions are shown as
solid black lines.
tain the desired degree corresponds to a column that
stretches over a row that has already been completed.
This situation is illustrated in Fig. 9. The routine
could start over until the adjacency matrix is success-
fully built, but such a strategy is inefficient. Instead,
after completion of each row, all uncompleted rows
are checked for possible future conflicts. If a conflict
is detected, entries in the row that causes the conflict
are shuffled randomly until the conflict is resolved.
If the conflict cannot be resolved after a predefined
number of iterations, the routine starts over. With
this strategy, more than 80% of the attempts are suc-
cessful, and the success rate becomes even higher for
larger matrix sizes, which we deem acceptable for our
purposes.
With this routine, bands of arbitrary shape can
be realised, as long as the bandwidth is larger than
the desired constant degree. To determine the node-
dependent bandwidth b according to Eqs. (12) and
(15), we determine ho
j
|H
†
OH − O|o
j
i for all eigen-
states |o
j
i of the observable and pick the maximum
value. Since this procedure becomes computationally
costly, we fit ∆o as an exponential function of L to
data for system sizes up to L = 10, and then ex-
trapolate to larger system sizes; see Fig. 10 for an
illustration.
Using b computed according to this description, to-
gether with the constant degree (10), defines the BRF
Figure 12: Linear fits of the standard deviations of the dis-
tributions of matrix elements H
jk
shown in Fig. 11.
ensemble. In the BRC ensemble, a constant band-
width max
j
b
o
(j) is used for all nodes. The BVF
ensemble also uses the exact bandwidth b, but con-
flicts as described above that may occur in the process
of creating the adjacency matrices are not resolved.
As a subsequent step after completing the last row,
excess edges originating from not resolved conflicts
are dropped randomly until the final matrix has the
same total number of nonzero elements as the original
Hamiltonian. As a result, the degree varies along the
nodes for graphs of that ensemble, but the average
degree agrees with the constant degree of the graphs
in the EXH ensemble.
Figure 13: Plots of adjacency matrices A as in Fig. 1, but
for magnetisation M with equal weights a
m
(i) ≡ 1. The
sorting of nodes according to the degenerate eigenvalues of
M is not unique.
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 14
Figure 14: As in Fig. 5, but using the magnetisation (9) with equal weights a
m
(i) ≡ 1 as observable.
D.2 Weight distributions of random graphs
The adjacency matrices obtained in Sec. D.1 charac-
terise the topology of the graph, but not the weights
of the edges. The random coupling coefficients of the
original Hamiltonian (8) cause the weights to be ran-
dom variables which we assume here to be indepen-
dent and identically distributed (i.i.d.). Since each
weight is a sum of many terms involving the coupling
coefficients, the central limit theorem can be expected
to apply, rendering the weights normally distributed.
We verify this expectation by numerically computing
the weights H
jk
. The histograms of weights shown
in Fig. 11 confirm Gaussian weight distributions with
zero mean. The standard deviation of the diagonal en-
tries, shown in Fig. 12, scale linearly with L, which en-
sures extensivity of the Hamiltonian, while the stan-
dard deviation of the off-diagonal entries remains con-
stant.
To impose the normalisation kHk = 1 onto the
various random graph ensembles we numerically de-
termine kHk for several system sizes between 4 and
10. Figure 10 shows that kHk increases linearly with
L, and we use a fit to these data for extrapolating
kHk to larger system sizes. For system sizes that are
beyond the reach of exact diagonalisation, this pro-
cedure is necessary because no algorithms are known
that can exploit the sparsity of H
jk
to efficiently de-
termine spectral properties such as the operator norm.
E Non-randomised magnetisation
In the present paper we mostly used as our observ-
able of choice a randomised version of the magnetisa-
tion M , with normally distributed random numbers
a
m
(i) in Eq. (9). This choice is motivated by the
convenience of working with an observable that has
a nondegenerate spectrum, such that the sorting of
eigenstates according to their eigenvalues is unique.
In this appendix we use the standard magnetisation
with equal weights a
m
(i) ≡ 1. We discuss the differ-
ence in the adjacency matrices brought about by this
choice of weights, and we demonstrate that the tools
and ideas put forward in this paper can nonetheless
be applied in this case.
In Fig. 13 we show the adjacency matrices for all
choices of n and d for a system of L = 8 spins. Unlike
in Fig. 1, the bandwidth is block shaped, where each
block corresponds to one of the 2L + 1 degenerate
eigenvalues of M, and the block size is given by the
degeneracy of the respective eigenvalue. The order of
the nodes corresponding to one eigenvalue is set by
the order the eigenstates are determined in python,
for which no pattern is apparent. Any node with the
eigenvalue closest to zero could work as an equilibrium
node.
To estimate the equilibration times, we calculate
the maximum flow value between the two nodes with
largest distance to equilibrium (k = 1 and k = N).
Since these eigenspaces are nondegenerate, the use of
the magnetisation with equal weights does not lead
to ambiguities. In Fig. 14 we show the equilibration
times T
eq
from the exact time evolution vs. the maxi-
mum flow value f
max
. The results are similar to those
for a randomised magnetisation in Fig. 5, showing
strong anticorrelations between t
eq
and f
max
.
References
[1] J. von Neumann. Beweis des Ergoden-
satzes und des H-Theorems in der neuen
Mechanik. Z. Phys., 57:30–70, 1929. DOI:
10.1007/BF01339852.
[2] H. Tasaki. From quantum dynamics to the
canonical distribution: General picture and a rig-
orous example. Phys. Rev. Lett., 80:1373–1376,
Feb 1998. DOI: 10.1103/PhysRevLett.80.1373.
[3] N. Linden, S. Popescu, A. J. Short, and A. Win-
ter. Quantum mechanical evolution towards ther-
mal equilibrium. Phys. Rev. E, 79:061103, 2009.
DOI: 10.1103/PhysRevE.79.061103.
[4] S. Goldstein, J. L. Lebowitz, C. Mastrodonato,
R. Tumulka, and N. Zanghì. Approach to ther-
mal equilibrium of macroscopic quantum sys-
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 15
tems. Phys. Rev. E, 81:011109, 2010. DOI:
10.1103/PhysRevE.81.011109.
[5] P. Reimann. Canonical thermalization. New
J. Phys., 12:055027, 2010. DOI: 10.1088/1367-
2630/12/5/055027.
[6] P. Reimann and M. Kastner. Equilibration of
isolated macroscopic quantum systems. New
J. Phys., 14:043020, 2012. DOI: 10.1088/1367-
2630/14/4/043020.
[7] C. Gogolin and J. Eisert. Equilibration, ther-
malisation, and the emergence of statistical me-
chanics in closed quantum systems. Rep. Prog.
Phys., 79:056001, 2016. DOI: 10.1088/0034-
4885/79/5/056001.
[8] T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda.
Thermalization and prethermalization in isolated
quantum systems: a theoretical overview. J.
Phys. B, 51:112001, 2018. DOI: 10.1088/1361-
6455/aabcdf.
[9] A. J. Short and T. C. Farrelly. Quantum equili-
bration in finite time. New J. Phys., 14:013063,
2012. DOI: 10.1088/1367-2630/14/1/013063.
[10] S. Goldstein, T. Hara, and H. Tasaki. Time scales
in the approach to equilibrium of macroscopic
quantum systems. Phys. Rev. Lett., 111:140401,
2013. DOI: 10.1103/PhysRevLett.111.140401.
[11] S. Goldstein, T. Hara, and H. Tasaki. Ex-
tremely quick thermalization in a macroscopic
quantum system for a typical nonequilibrium
subspace. New J. Phys., 17:045002, 2015. DOI:
10.1088/1367-2630/17/4/045002.
[12] A. S. L. Malabarba, L. P. García-Pintos, N. Lin-
den, T. C. Farrelly, and A. J. Short. Quantum
systems equilibrate rapidly for most observables.
Phys. Rev. E, 90:012121. DOI: 10.1103/Phys-
RevE.90.012121.
[13] T. Farrelly. Equilibration of quantum gases. New
J. Phys., 18:073014, 2016. DOI: 10.1088/1367-
2630/18/7/073014.
[14] P. Reimann. Typical fast thermaliza-
tion processes in closed many-body sys-
tems. Nat. Commun., 7:10821, 2016. DOI:
10.1038/ncomms10821.
[15] H. Wilming, M. Goihl, C. Krumnow, and J. Eis-
ert. Towards local equilibration in closed in-
teracting quantum many-body systems. URL
https://arxiv.org/abs/1704.06291.
[16] L. P. García-Pintos, N. Linden, A. S. L. Mal-
abarba, A. J. Short, and A. Winter. Equi-
libration time scales of physically relevant ob-
servables. Phys. Rev. X, 7:031027, 2017. DOI:
10.1103/PhysRevX.7.031027.
[17] T. R. de Oliveira, C. Charalambous,
D. Jonathan, M. Lewenstein, and A. Riera.
Equilibration time scales in closed many-body
quantum systems. New J. Phys., 20:033032,
2018. DOI: 10.1088/1367-2630/aab03b.
[18] P. Reimann. Transportless equilibration in iso-
lated many-body quantum systems. New J.
Phys., 21:053014, 2019. DOI: 10.1088/1367-
2630/ab1a63.
[19] E. P. Wigner. Characteristic vectors of bordered
matrices with infinite dimensions. Ann. Math.,
62:548–564, 1955. DOI: 10.2307/1970079.
[20] E. P. Wigner. Characteristic vectors of bordered
matrices with infinite dimensions II. Ann. Math.,
65:203–207, 1957. DOI: 10.2307/1969956.
[21] E. P. Wigner. On the distribution of the roots
of certain symmetric matrices. Ann. Math., 67:
325–327, 1958. DOI: 10.2307/1970008.
[22] M. L. Mehta. Random Matrices. Elsevier, Ams-
terdam, 3rd edition, 2004.
[23] L. F. Santos and E. J. Torres-Herrera. Ana-
lytical expressions for the evolution of many-
body quantum systems quenched far from equi-
librium. AIP Conf. Proc., 1912:020015, 2017.
DOI: 10.1063/1.5016140.
[24] E. J. Torres-Herrera, A. M. García-García, and
L. F. Santos. Generic dynamical features of
quenched interacting quantum systems: Survival
probability, density imbalance, and out-of-time-
ordered correlator. Phys. Rev. B, 97:060303,
2018. DOI: 10.1103/PhysRevB.97.060303.
[25] L. F. Santos and E. J. Torres-Herrera. Nonequi-
librium many-body quantum dynamics: From
full random matrices to real systems. In
F. Binder, L. Correa, C. Gogolin, J. Anders,
and G. Adesso, editors, Thermodynamics in
the Quantum Regime. Fundamental Theories of
Physics, pages 457–479. Springer, 2018. DOI:
10.1007/978-3-319-99046-0_19.
[26] E. J. Torres-Herrera, J. Karp, M. Távora, and
L. F. Santos. Realistic many-body quantum sys-
tems vs. full random matrices: Static and dy-
namical properties. Entropy, 18:1–20, 2016. DOI:
10.3390/e18100359.
[27] J. B. French and S. S. M. Wong. Validity of
random matrix theories for many-particle sys-
tems. Phys. Lett. B, 33:449–452, 1970. DOI:
10.1016/0370-2693(70)90213-3.
[28] O. Bohigas and J. Flores. Two-body ran-
dom Hamiltonian and level density. Phys.
Lett. B, 34:261–263, 1971. DOI: 10.1016/0370-
2693(71)90598-3.
[29] V. K. B. Kota, A. Relaño, J. Retamosa, and
M. Vyas. Thermalization in the two-body ran-
dom ensemble. J. Stat. Mech., 2011:P10028,
2011. DOI: 10.1088/1742-5468/2011/10/P10028.
[30] V. V. Flambaum and F. M. Izrailev. Entropy
production and wave packet dynamics in the
fock space of closed chaotic many-body sys-
tems. Phys. Rev. E, 64:036220, 2001. DOI:
10.1103/PhysRevE.64.036220.
[31] F. Borgonovi, F. M. Izrailev, and L. F. Santos.
Exponentially fast dynamics of chaotic many-
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 16
body systems. Phys. Rev. E, 99:010101, 2019.
DOI: 10.1103/PhysRevE.99.010101.
[32] F. Borgonovi and F. M. Izrailev. Emergence of
correlations in the process of thermalization of in-
teracting bosons. Phys. Rev. E, 99:012115, 2019.
DOI: 10.1103/PhysRevE.99.012115.
[33] V. K. B. Kota. Embedded random matrix ensem-
bles for complexity and chaos in finite interacting
particle systems. Phys. Rep., 347:223–288, 2001.
DOI: 10.1016/S0370-1573(00)00113-7.
[34] V. K. B. Kota and N. D. Chavda. Random k-
body ensembles for chaos and thermalization in
isolated systems. Entropy, 20:1–22, 2018. DOI:
10.3390/e20070541.
[35] M. Vyas and T. H. Seligman. Random ma-
trix ensembles for many-body quantum systems.
AIP Conf. Proc., 1950:030009, 2018. DOI:
10.1063/1.5031701.
[36] T. H. Seligman, J. J.M. Verbaarschot, and M. R.
Zirnbauer. Spectral fluctuation properties of
Hamiltonian systems: The transition region be-
tween order and chaos. J. Phys. A, 18:2751–2770,
1985. DOI: 10.1088/0305-4470/18/14/026.
[37] Y. V. Fyodorov, O. A. Chubykalo, F. M. Izrailev,
and G. Casati. Wigner random banded matri-
ces with sparse structure: Local spectral density
of states. Phys. Rev. Lett., 76:1603–1606, 1996.
DOI: 10.1103/PhysRevLett.76.1603.
[38] L. Erdos and A. Knowles. Quantum diffusion and
eigenfunction delocalization in a random band
matrix model. Commun. Math. Phys., 303:509–
554, 2011. DOI: 10.1007/s00220-011-1204-2.
[39] A. D. Mirlin, Y. V. Fyodorov, F.-M. Dittes,
J. Quezada, and T. H. Seligman. Transition from
localized to extended eigenstates in the ensem-
ble of power-law random banded matrices. Phys.
Rev. E, 54:3221–3230, 1996. DOI: 10.1103/Phys-
RevE.54.3221.
[40] P. A. Nosov and I. M. Khaymovich. Robust-
ness of delocalization to the inclusion of soft
constraints in long-range random models. Phys.
Rev. B, 99:224208, 2019. DOI: 10.1103/Phys-
RevB.99.224208.
[41] J. A. Méndez-Bermúdez, G. F. De Arruda, F. A.
Rodrigues, and Y. Moreno. Diluted banded ran-
dom matrices: Scaling behavior of eigenfunction
and spectral properties. J. Phys. A, 50:495205,
2017. DOI: 10.1088/1751-8121/aa9509.
[42] I. Jana and A. Soshnikov. Distribution of singu-
lar values of random band matrices; Marchenko–
Pastur law and more. J. Stat. Phys., 168:964–
985, 2017. DOI: 10.1007/s10955-017-1844-5.
[43] P. Bourgade. Random band matrices. Proc.
Int. Cong. Math., 3:2745–2770, 2018. DOI:
10.1142/11060.
[44] I. Dumitriu and Y. Zhu. Sparse general Wigner-
type matrices: Local law and eigenvector delocal-
ization. J. Math. Phys., 60:023301, 2019. DOI:
10.1063/1.5053613.
[45] F. Borgonovi, F. M. Izrailev, L. F. Santos,
and V. G. Zelevinsky. Quantum chaos and
thermalization in isolated systems of interact-
ing particles. Phys. Rep., 626:1–58, 2016. DOI:
10.1016/j.physrep.2016.02.005.
[46] G. P. Brandino, A. De Luca, R. M. Konik, and
G. Mussardo. Quench dynamics in randomly
generated extended quantum models. Phys.
Rev. B, 85:214435, 2012. DOI: 10.1103/Phys-
RevB.85.214435.
[47] S. Goldstein, J. L. Lebowitz, R. Tumulka, and
N. Zanghì. Long-time behavior of macroscopic
quantum systems. Eur. Phys. J. H, 35:173–200,
2010. DOI: 10.1140/epjh/e2010-00007-7.
[48] D. Nickelsen and M. Kastner. Classical
Lieb-Robinson bound for estimating equilibra-
tion timescales of isolated quantum systems.
Phys. Rev. Lett., 122:180602, 2019. DOI:
10.1103/PhysRevLett.122.180602.
[49] R. Diestel. Graph Theory. Springer, Berlin, 3rd
edition, 2005.
[50] I. Arad, T. Kuwahara, and Z. Landau. Con-
necting global and local energy distributions in
quantum spin models on a lattice. J. Stat.
Mech., 2016:033301, 2016. DOI: 10.1088/1742-
5468/2016/03/033301.
[51] W. Beugeling, R. Moessner, and M. Haque.
Off-diagonal matrix elements of local opera-
tors in many-body quantum systems. Phys.
Rev. E, 91:012144, 2015. DOI: 10.1103/Phys-
RevE.91.012144.
[52] J. Wang and W.-G. Wang. Correlations in eigen-
functions of quantum chaotic systems with sparse
Hamiltonian matrices. Phys. Rev. E, 96:052221,
2017. DOI: 10.1103/PhysRevE.96.052221.
[53] S. C. Morampudi and C. R. Laumann. Many-
body systems with random spatially local inter-
actions. Phys. Rev. B, 100:245152, 2019. DOI:
10.1103/PhysRevB.100.245152.
[54] C. Nation and D. Porras. Off-diagonal observ-
able elements from random matrix theory: dis-
tributions, fluctuations, and eigenstate thermal-
ization. New J. Phys., 20:103003, 2018. DOI:
10.1088/1367-2630/aae28f.
[55] L. Foini and J. Kurchan. Eigenstate thermal-
ization hypothesis and out of time order cor-
relators. Phys. Rev. E, 99:42139, 2019. DOI:
10.1103/PhysRevE.99.042139.
[56] M. Brenes, S. Pappalardi, J. Goold, and A. Silva.
Multipartite entanglement structure in the eigen-
state thermalization hypothesis. Phys. Rev.
Lett., 124:040605, 2020. DOI: 10.1103/Phys-
RevLett.124.040605.
Accepted in Quantum 2020-05-21, click title to verify. Published under CC-BY 4.0. 17