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 [18]. 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 [918]. 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 [1921]. 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 [2325]. 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 [2932]. General-
isations of the two-body random ensemble to random
k-body interactions have been studied for fermionic
as well as bosonic systems [3335]. 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 [3638]. 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 [4244]. 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) =
n1
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
±
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
jk,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 = (2qL)/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 [5153], a topic that
has been discussed also as an extension to the
eigenstate thermalisation hypothesis [5456]. 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
18L
. 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 pair