Modelling equilibration of local many-body quantum systems by random graph ensembles

We introduce structured random matrix ensembles, constructed to model many-body quantum systems with local interactions. These ensembles are employed to study equilibration of isolated many-body quantum systems, showing that rather complex matrix structures, well beyond Wigner's full or banded random matrices, are required to faithfully 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 network theory. One of these tools, called the maximum flow value, is found to be an excellent proxy for equilibration times. Since maximum flow values are less expensive to compute, they give access to approximate equilibration times for system sizes beyond those accessible by exact diagonalisation.


Introduction
By means of typicality arguments, equilibration of isolated quantum systems has been rigorously proved to occur under rather general conditions [1][2][3][4][5][6][7][8]. It remains however a key open problem to identify which properties of quantum systems are responsible for ensuring that equilibration timescales are physically realistic [9][10][11][12][13][14][15][16][17][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 certain properties of the nuclei [19][20][21]. One such choice, the so-called Gaussian ensembles, is given by matrices where all elements are normally distributed i.i.d. numbers normalized to some measure [22]. Gaussian random matrices are a particularly convenient choice, as they allow for the analytic calculation of ensemble-averaged properties like limiting distributions 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][24][25]. Since the time evolution of expectation values of observables is governed by differences of energy eigenvalues, the level-spacing statistics of random matrices influences the equilibration times [26]. The suppression of small spacings, as in the Wigner-Dyson distribution, favours short equilibration 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 suggests that Gaussian random matrices are unsuitable for representing physically realistic quantum manybody Hamiltonians.
Locality, in the sense of short-ranged interactions 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 fullyconnected systems with distance-independent interaction 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 random 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 random ensemble has also been used to study thermalisation of isolated fermionic systems [29][30][31][32]. Generalisations of the two-body random ensemble to random k-body interactions have been studied for fermionic as well as bosonic systems [33][34][35]. More generally, random matrix ensembles that use full random matrices 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 matrices and that they display a banded structure. Similarly, the (non-random) Hamiltonian of a manybody system with local interactions, when represented in the tensor product basis with respect to which locality is defined, turns out to be a banded, sparse matrix. In fact, already Wigner in his early attempts to model nuclei of heavy atoms introduced banded matrices that have nonzero elements only up to a certain distance from the diagonal [19,20]. In addition to resembling Hamiltonians with local interactions, several studies have shown that random band matrices are also capable of modelling certain aspects important for equilibration timescales, e.g. the transition from Wigner-Dyson statistics for chaotic system to Poisson statistics of regular systems, and the emergence of many-body localisation [36][37][38]. Variants of Wigner band matrices include power-law banded random 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 physical systems more realistically, they come with the disadvantage that analytic calculations of their spectral properties are far more intricate [42][43][44]. An overview of applications of embedded and banded ensembles is provided in Ref. [45]. In a different attempt to include aspects of locality into random Hamiltonians, the authors of Ref. [46] disregarded the banded structure 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 (bandedness, sparseness, etc.) of its matrix representation. To this purpose, we interpret the Hamiltonian matrix as the connectivity matrix of a weighted graph, which allows us to apply concepts from graph and network theory. As our first main result, we identify key features of the graphs corresponding to local Hamiltonians, including degree distributions as well as a node-dependent ("curved") bandedness of the connectivity 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 either the correct bandshape or constant degree, however, 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 calculate its maximum flow value, from which the equilibration time can be deduced. Combining the efficient graph construction and the calculation of the maximum flow value, equilibration times of local Hamiltonians can be investigated more efficiently.

Equilibration
We choose the equilibration dynamics of a quantum 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 generated by a time-independent Hamiltonian H. Under the dynamics of (1), the infinite time average of the expectation value of an observable O can be written as where |e j are energy eigenstates and c j = e j |ψ 0 are the overlaps with the initial state. Unlike the microcanonical ensemble that uses equal weights, the socalled 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 O 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 O evolves quasi-periodically in time. As a consequence, even at arbitrarily late times, fluctuations around the equilibrium value O diag occur. Because of these fluctuations, a sensible definition of equilibration in isolated macroscopic quantum systems requires a probabilistic perspective. The largest part of the Hilbert space, referred to as the equilibrium subspace, is spanned by eigenstates |o j of O that have eigenvalues o j close to the equilibrium value, o j ≈ O diag . The complement of this subspace is the nonequilibrium subspace, spanned by eigenstates of O with eigenvalues that appreciably differ from O 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 |ψ not only guarantees that ψ|O|ψ ≈ O 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 according to the Haar measure will most likely be an equilibrium 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 becomes vanishingly small with increasing system size.
Throughout this section, equilibration was discussed with reference to a specific observable O, and also the (non)equilibrium subspaces have been defined 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 finitedimensional 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 e j | onto energy eigenspaces), or that evolve periodically and hence do not equilibrate (like |e i e j | + |e j e 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 respect to specific observables is hence not a shortcoming, but a desired feature.

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 of O may play a distinguished role. In this basis, the Schrödinger equation (1) readṡ where The eigenstates |o j are assumed to be sorted in ascending order, o j+1 ≥ o j . By componentwise integration, Eq. (3) can be written as 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 equilibration 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 amplitudes x j (t), oscillations need to propagate through the network in order to activate the equilibrium oscillators. Once these oscillators are active, measurement probabilities x j for equilibrium values o j ≈ O 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 distribution, 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 matrix A with matrix elements We use 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) ≈ 0 are nonequilibrium nodes. In Ref. [48], this network picture of quantum dynamics has been used to derive a lower bound on the timescale at which the system equilibrates. Remarkably, the bound was found to depend on the degrees of locality of the Hamiltonian H as well as of the observable 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 corresponding 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 equilibration 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 oscillation of a nonequilibrium node propagating through the network to an equilibrium node along all of the existing pathways between the nodes, promising candidates are measures that take into account multiple pathways. The node connectivity, defined as the minimum number of nodes that need to be removed in order to split the graph into two disconnected subgraphs, 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 network, where f max is the total capacity between two nodes. To define f max [49], we denote the flow from node u to node v with f (u, v) and impose the capacity constraint f (u, v) < |H uv | for all u and v. The flow value is defined as |f (u)| = v f (u, v), that is the sum of all flows leaving node u. The maximum flow value f max (k, j) is then the maximum of |f (k)| under the additional constraint that |f (k)| = |f (j)|, in other words, imposing that all flow entering the network at k must be able to reach j. 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. We will apply the maximum flow measure to the network picture of quantum spin models in Sec. 5.3.

Spin model
There exist various notions of locality, and our goal is to explore network features and exploit graph measures 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., statistical 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 interaction 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-dimensional 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 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 components 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 guarantees that accidental cancellations in the matrix elements H jk occur with zero probability, and hence the elements of the adjacency matrix (6) do not vanish 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 H = 1 to get rid of overall factors in H, as a Hamiltonian 2H would equilibrate twice as fast as H for trivial reasons. Unnormalised Hamiltonians usually scale extensively or superextensively with the lattice size L, and our normalisation convention will therefore lead to longer equilibration times compared to the unnormalised counterpart. We will comment on implications 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, with real weights a m . The eigenstates |m j 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 eigenvalues occur with high degeneracies, and hence the eigenstates |m 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 deviation 1/10. A unique order of eigenstates brings out the graph properties more clearly and demonstrates that our analysis works for generic observables. Lifting the degeneracy, however, is not a requirement for our analysis, as is demonstrated in Appendix E. Due to the independent, anisotropic and inhomogeneous coupling, we expect the equilibrium value M diag to be close to zero, which is confirmed by exact diagonalisation. Hence, equilibrium eigenstates are at the centre and nonequilibrium eigenstates at the edges of the spectrum.

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 magnetisation M (9). We focus in particular on the dependence of the graph properties on the locality parameters n and d introduced in (8).
The degree distribution is a prominent graph property. The degree (j) = k =j A jk of a node j counts the number of edges of this node. The degree distribution 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 Poisson 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 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 Gaussian orthogonal ensemble. Other choices of n and d create matrices of intermediate sparsity.
The plots in Fig. 1 reveal that, in addition to being sparse, the adjacency matrices of sufficiently local 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 matrices may be more suitable to model physical systems [19,20]. Also, a banded structure is consistent with 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 suppressed nonzero elements outside the band. The reasoning 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 exponentially small matrix elements outside the bands [48]. However, for the choice of Hamiltonian and observable 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 observe 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 diagonal 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 implies that a nonequilibrium node j is connected only to nodes with a similar index i ≈ j, whereas equilibrium nodes may also be connected to nodes whose index i differs more substantially from j. For the adjacency matrix of the spin Hamiltonian (8), the nodedependent bandwidth can be calculated analytically. We first derive such a result for the magnetisation (9) with homogeneous weights a m (i) ≡ 1, and then generalise the argument to arbitrary observables.
Writing the Hamiltonian in terms of ladder operators σ ± i = 1 2 (σ x i ± iσ y i ), n-spin terms translate into products of n ladder operators, which allow for at most n spin flips and can change the magnetisation by at most ∆m = 2n/L. We define the bandwidth 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,    (12) and applies to the case am(i) ≡ 1 for the magnetisation M in (9) as observable. The bandwidth bm(j) follows from the general formula (14), here applied to the randomised M using normaldistributed am(i) as defined in (9). The bandwidthsb(j) andbm(j) follow from the definition in (13) for non-randomised and randomised M respectively.
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, 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, Here, g(o) is the density of eigenstates O, and G(o) the integrated density of states, defined as the number of eigenstates 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 of O. Equation (14), being general, also applies to the magnetisation (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 actual structure of the adjacency matrix A jk of the spin model (8) when using the magnetisation (9) as observable. Overall, we find that the locality parameter n determines the bandwidth of the adjacency matrix, whereas n and d determine the density of nonzero elements inside the band.

Comparison of random graph ensembles
In the previous section we have studied certain graph features of adjacency matrices, in particular their degree distributions and node-dependent bandwidths b, and investigated the dependence of these features on the locality parameters n and d of the Hamiltonian. Now we reverse the logic by constructing random graphs with a certain b and , and investigate to what extend the corresponding random matrix ensembles reproduce the equilibration dynamics of Hamiltonians with locality parameters n and d. In the following we sketch the main features of the random graph ensembles considered in this paper. BRC: Banded regular graphs with constant (nodeindependent) bandwidth b ≡ max j b m (j) using (14) are generated randomly but not uniformly. The maximum bandwidth and the constant degree 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 ensembles have the same total number of edges, i.e. all resulting Hamiltonians have the same number of nonzero matrix elements. Apart from the EXH ensemble, the weights of all graph edges are i.i.d. numbers sampled from normal distributions matching the statistics of the exact H jk . Details of the construction of these random graphs are given in Appendix D.
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 Hamiltonians and assess which of the ensembles lead to timescales that faithfully reproduce those of Hamiltonians 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 H = 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 magnetisation.
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 O and O 2 for the first time reach their diagonal averages with an allowed margin of error of 10% of their respective initial 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 approximately 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) increase 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 oscillations 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, differ, 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 between the random matrix elements of the Hamiltonian and the observable [51][52][53], a topic that has been discussed also as an extension to the eigenstate thermalisation hypothesis [54][55][56]. To fully focus on the effects of graph topology, it may make sense to compare equilibration times 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 best suited to reproduce the equilibration timescales of the ensembles EXH and EXA. However, on the basis of the available data, one may speculate that the agreement deteriorates with decreasing locality, as in the right plot of Fig. 4.
We conclude that constant degree and bandedness are essential features of the adjacency matrices of local Hamiltonians that need to be taken into account in a random matrix approach. The shape of the nodedependent bandwidth also has a role to play, but appears to be somewhat less crucial. Since it is computationally costly to calculate the exact adjacency matrix required for the EXA ensemble, resorting to BRF as the next-best, and computationally less expensive, ensemble is a potential avenue for studying equilibration times of larger systems in a random matrix setting. 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 constant (system-size independent) timescale: because of the homogeneity, no currents across macroscopic distances 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 H = 1, however, introduces an extra factor of 1/ H in the time-evolution operator, resulting in the observed approximately linear increase of T eq with L.

Maximum flow vs. equilibration time
To determine the equilibration time for a given Hamiltonian, 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 calculating 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 systems, we fix n and d and calculate f max and T eq for random samples of Hamiltonians for system sizes between 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 linear fit to the data, shown in the plot, is used in Sec. 5.4 to estimate equilibration times of larger systems.

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

Conclusions
We introduced and studied various random graph ensembles, with the aim of modelling equilibration in a spin system with a given degree of locality. Interpreting 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 rational of our analysis proceeds in three steps: (a) Working with a local Hamilton H and a local observable O implies the existence of a basis in which H and O are simultaneously sparse. A numerical study of quantum spin models with random 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 numerically study samples of Hamiltonians from the various ensembles, finding that the BRF ensemble faithfully captures the dependence of equilibration times on locality, while having computational 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 correlated with the equilibration time T eq of the corresponding 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 integration.
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 single 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 equilibration times for large system sizes can then be studied using f max of the generated graphs, without deriving 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 equilibration 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 described 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 connectivity and the shortest path length, but these turned out to be less strongly correlated with T eq . Finally, we believe that the random graph ensembles 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.
A More on the network perspective In Ref. [48] it has been argued that the quantity 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 its maximal values, we can say the system has equilibrated.
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 (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.
and hence Eq. 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 = (−i H t) q /q!. For t = 0, nodes j and k are uncoupled. For small times t < 1/ H , |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/ H , 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 statistics 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 sufficiently large t. This is the saturation mentioned above, which indicates that the system has equilibrated.
The inverse of the norm of the Hamiltonian, 1/ H , therefore gives an indication on the equilibration timescale, but it does not capture the influence of locality. Locality enters through the vanishing coupling chains of small lengths q. For the dense Hamiltonian of a nonlocal system, the weights w q on short coupling lengths q ≈ 2 for small t are already sufficient 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 nearestneighbour 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 negligible for t ≈ H . 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 ≈ H .

B Constant degree
As a consequence of statistical isotropy and homogeneity of the Hamiltonian, all nodes of a given adjacency graph have a constant degree . By numerically determining for system sizes up to L = 16 and various combinations of locality parameters n and d, we conjectured that where the first term depends linearly on L, and the second, not yet determined second term C(n, d), is independent of L. By studying C for various fixed values of n we were able to further conjecture that In the last step of this calculation we used that is the number of distinct ways to draw a times from b options with replacement where the order 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 which agrees with Eq. (10). We tested (10) for all allowed combinations of locality parameters n and d for systems 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 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 coordinate 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 function 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.
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 efficiently 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 degree . 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 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 bottom, 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. from all possible matrices. We are not aware of similar 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 corresponding 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 number e of nonzero elements already placed in each uncompleted row, and choose only as many nonzero entries 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, reducing 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 completed 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 obtain 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 successfully 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 successful, 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 nodedependent bandwidth b according to Eqs. (12) and (15), we determine o j |H † OH − O|o j for all eigenstates |o j 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 extrapolate to larger system sizes; see Fig. 10 for an illustration.
Using b computed according to this description, together with the constant degree (10), defines the BRF ensemble. In the BRC ensemble, a constant bandwidth max j b o (j) is used for all nodes. The BVF ensemble also uses the exact bandwidth b, but conflicts 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.

D.2 Weight distributions of random graphs
The adjacency matrices obtained in Sec. D.1 characterise 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 random variables which we assume here to be independent 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 entries, shown in Fig. 12, scale linearly with L, which ensures extensivity of the Hamiltonian, while the standard deviation of the off-diagonal entries remains constant.
To impose the normalisation H = 1 onto the various random graph ensembles we numerically determine H for several system sizes between 4 and 10. Figure 10 shows that H increases linearly with L, and we use a fit to these data for extrapolating H to larger system sizes. For system sizes that are beyond the reach of exact diagonalisation, this procedure is necessary because no algorithms are known that can exploit the sparsity of H jk to efficiently determine spectral properties such as the operator norm.

E Non-randomised magnetisation
In the present paper we mostly used as our observable of choice a randomised version of the magnetisation 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 difference 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 maximum 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 .