The Thermodynamic Limit of Spin Systems on Random Graphs

We utilise the graphon--a continuous mathematical object which represents the limit of convergent sequences of dense graphs--to formulate a general, continuous description of quantum spin systems in thermal equilibrium when the average co-ordination number grows extensively in the system size. Specifically, we derive a closed set of coupled non-linear Fredholm integral equations which govern the properties of the system. The graphon forms the kernel of these equations and their solution yields exact expressions for the macroscopic observables in the system in the thermodynamic limit. We analyse these equations for both quantum and classical spin systems, recovering known results and providing novel analytical solutions for a range of more complex cases. We supplement this with controlled, finite-size numerical calculations using Monte-Carlo and Tensor Network methods, showing their convergence towards our analytical results with increasing system size.

Introduction -The physical properties of interacting systems are strongly affected by the connectivity of their components.For instance, network topology plays a decisive role in the rate of disease spreading in infectious disease models [1] whilst systematic studies have been undertaken into the affect of connectivity on the synchronisation of oscillators [2][3][4][5].
In interacting spin systems, the same ideas hold true: frustration causes the manifestation of exotic phases of matter such as a spin liquids [6] and small-world effects alter the underlying universality class of ordereddisordered phase transitions [7].The difficulty of solving the many-body problem (especially in the quantum regime), however, means a more general characterisation of how network topology influences strongly correlated systems is unknown.
When disorder is absent and the average co-ordination number becomes large, interacting many-body systems fall into the mean-field universality class and become amenable to simpler, mathematical and computational approaches [8][9][10][11][12].Despite mostly being applied to translationally-invariant systems, the mean field approach is known to be valid for an infinite multitude of networks, whether homogeneous or heterogeneous [13].Whilst it is only exact in the thermodynamic limit and when the average co-ordination number grows proportionally with the system size [14], mean-field theory can provide meaningful physical predictions for lowdimensional systems [15,16].Within the field of graph theory, network structures-whether heterogeneous or homogeneous-with an extensive co-ordination number are well-characterised.Their thermodynamic limit is succinctly described by the graphon [17,18], a continuous mathematical object which represents the limit of a sequence of adjacency matrices as the number of vertices tends to infinity and the average co-ordination number grows extensively.* jtindall@flatironinstitute.org Here we utilise the graphon in the study of interacting spin systems.This allows us to formally take the thermodynamic limit and derive an exact, continuous theory for the limit of sequences of discrete Hamiltonian on graphs of increasing size and average co-ordination number.Specifically, we take a very general spin Hamiltonian defined over an arbitrary graph and, for sequences of dense graphs whose limit is known to converge to a given graphon, derive a coupled set of integral equations which exactly describe the equilibrium physics of the limit of the corresponding sequence of Hamiltonians.The graphon forms the kernel in these integral equations and the Physics of the system can be directly studied as a function of this object.
Taking several classical and quantum example models we demonstrate the utility of these integral equations: i) verifying previous results on all-to-all spin systems, ii) proving the existence of a finite-temperature phase transition in the classical Ising model for any graphon and iii) deriving analytical solutions for the equilibrium observables of spin models on novel non-trivial, heterogeneous networks.We reinforce our analytical solutions with large-scale, finite-size Monte Carlo and Tensor Network simulations.Whilst the spin systems we treat in this work are commonly studied due to their relevance as models of real-world magnetism, they also find application in many other branches of science, including the political, social and biological sciences [19][20][21].
Hamiltonian -Our starting point is L qubits placed on the L vertices of a graph G L .The graph is specified by an L×L symmetric adjacency matrix A G L with elements A v,v ′ which dictate the (weighted) connections between vertices (qubits) v, v ′ ∈ [1...L], [1..

.L]. The Hamiltonian reads
along the α-spin-axis on vertex v. Our focus is on graphs where 0 ) then we refer to the graph as 'dense' because the the average co-ordination number diverges with system size.The factor of 1 L in H(G L ) is necessary to ensure a finite, non-trivial energy density.A variety of well-known models-including the Curie-Weiss [22] and Lipkin Meshkov Glick [23,24] models-are contained within our Hamiltonian.The restriction 0 ≤ A v,v ′ ≤ 1 and the homogeneous nature of the field strengths, however, precludes our Hamiltonian from including disordered models such as spin-glass systems [25,26].
Theory -In this work, by utilising tools from graph theory and mean-field theory, we formulate an explicit, exact, continuous description of this system in thermal equilibrium in the thermodynamic limit.In order to describe our continuum formalism we must first introduce the concept of a graphon.This can be done by taking the vertices v = 1...L of a graph G L and performing the change of variables: , a real symmetric stepped function over the unit square such that for a given (x, y) ∈ I v × I v ′ , where Equipped with a well-defined metric for the similarity of two graphs, it can be shown that for certain sequences of graphs (G L ) L∈N the limit lim L→∞ W G L (x, y) converges to a well-defined symmetric function W (x, y) known as the 'graphon' [17,18].In the Supplemental Material (SM) we discuss these metrics in detail and provide theorems on the convergence of graph sequences under these metrics.Importantly, it is also possible to move in the opposite direction and given a graphon W (x, y) construct sequences of finite graphs whose limit is W (x, y).These finite graphs can be constructed via one of two methods: 'stochastic' or 'weighted' sampling of W (x, y) and we use G S L and G W L to refer to their respective realisation over the vertices v = 1....L.They can be constructed by defining the quantity (2) The adjacency matrix of the unweighted graph G S L is then defined by setting A v,v ′ = 1 with probability P v,v ′ and A v,v ′ = 0 otherwise.The adjacency matrix of the weighted graph G W L is defined by setting A v,v ′ = P v,v ′ .A given sequence of such realisations is guaranteed to converge to the graphon W (x, y) in the limit L → ∞ [27].
With the definition of the graphon in hand, the central result of this paper can be presented.
..) be a sequence of finite-size graphs generated as stochastic or weighted realisations of the graphon W (x, y).Then for a given inverse temperature T = 1/β the macroscopic properties of the equilibrium states of the sequence of Hamiltonians H(G L ) L∈N = (H(G 1 ), H(G 2 ), ...) converges and are determined by the solution of the following coupled integral equations 2 and the three functions λ α (x), with α ∈ {x, y, z}, each being continuous, real-valued and defined over the domain [0, 1].
In order to prove this theorem and arrive at Eq. ( 3) we state the following intermediate theorem L and G W L be the stochastic and weighted realisations on L vertices of a graphon W (x, y) respectively.For an arbitrary set of real, finite values for the parameters {J x , J y , J z , h x , h y , h z } the following is true This theorem is a significant generalisation of theorem 1 in Ref. [28] which proved this result solely for sequences of Erdős-Rényi graphs, which correspond to the constant graphon.The proof of theorem 2 (which can be found in the SM) relies on more general statistical properties of random graphs.With theorem 2 in hand, theorem 1 follows by: i) focusing strictly on the sequence (G W L ) L∈N = (G W 1 , G W 2 , ...) of weighted finite realisations of W (x, y), ii) applying meanfield theory (which is exact here in the thermodynamic limit) and iii) taking the continuum limit of the resulting equations by invoking the definition of the graphon.The SM contains full proofs of both theorem 1 and theorem 2.
If we can solve Eq.
(3) for the functions {λ x (x), λ y (x), λ z (x)}, then we have determined the equilibrium physics of the limit of the sequence (H(G 1 ), H(G 2 ), ...).The functions {λ x (x), λ y (x), λ z (x)} are a change of variables from the continuum limit of the spin degrees of freedom in the Hamiltonian.They directly encode the physical properties of the equilibrium state: the magnetisation on site v in the thermodynamic limit is specified by ⟨σ α (x)⟩ with x = lim L→∞ v L and is related to the λ functions by The total magnetisation along a given spin direction is The validity of the mean-field approximation here means we can compute multi-point correlators as products of on-site expectation values.How can we solve Eq.
In general there is no analytical solution and we will be restricted to numerical methods.Nonetheless, there are certain cases where they can be solved analytically.Consider the case the graphon is degenerate, i.e.W (x, y) where n is finite and f i (x) : [0, 1] → [0, 1].Substitution into the above equation tells us λ α (x) = n i=1 c α i f i (x) where c α i are real-valued coefficients which depend on the field strengths h α , couplings J α and the inverse temperature β but do not depend on x.These coefficients c α i are the solution of the set of 3n coupled equations which result from the substitution of λ α (x) = n i=1 c α i f i (x) into Eq.( 3).For a given set of J α , h α and value of β we therefore have a closed form for λ α (x) and various observables in the system.In our examples in the main text (further examples, including non-degenerate graphons are considered in the SM) we focus on n = 1 as they can be manipulated to yield closed forms for the equilibrium properties of the system.
Classical Ising model -We first set J x = J y = h x = h y = h z = 0 and J z = −1, realising the classical Ising model with zero field.Utilising sgn(z) tanh(β|z|) = tanh(z), our integral equations reduce to The Z 2 spin-flip symmetry is encoded in the fact that if λ z (x) is a solution to the equation then so is −λ z (x).Moreover, there is clearly always the trivial solution λ z (x) = 0 ∀x which corresponds to the disordered paramagnetic state with 0 magnetisation.Applying Banach's fixed-point theorem [29] to Eq. ( 7) tells us that, with certainty, when β < sup x∈[0,1] 1 0 W (x, y)dy this is the only solution.For larger values of β, however, there exists a non-trivial solution which corresponds to a ferromagnetic phase.For instance, when β → ∞ we have λ z (x) = 1 0 W (x, y)dx ̸ = 0 ∀x.Thus, following this analysis, we know that λ z (x, β) cannot be smooth and continuous over x ∈ [0, 1] and β ∈ [0, ∞] and there must exist a finite-order transition between the ferromagnetic solution and the paramagnetic solution at some critical temperature.Our continuum description has therefore allowed us to prove the existence of a ferromagnetic-paramagnetic phase transition for the Ising model on any dense graphwith a corresponding analytical upper bound on this temperature.A similar argument can be applied to a number of the limits of Eq. (1).Now let us treat some explicit examples.We first consider W (x, y) = p whose stochastic realisations are G ER (p): the Erdős-Rényi graph over L vertices where each edge appears independently with probability p. Observe from Eq. ( 7) that in this case λ z (x) = λ z = p tanh(βλ z ) and is independent of x.Substituting this into Eq.( 6) gives us the familiar self-consistent equation M z = tanh(βpM z ) for the magnetisation M z of the classical Ising model under the mean-field approximation.The edge probability p re-scales the temperature in the all-to-all model and the randomness of the model has no effect on the macroscopic physics in the thermodynamic limit-a result which has been proven to be general for spin systems on Erdős-Rényi graphs [22,28].
We consider the, more complex, separable graphon W (x, y) = xy, whose stochastic relatisations dictates that each pair of spins v and v ′ interacts with a strength 1 with probability (vv ′ /L 2 ) and strength 0 otherwise.One can also choose to directly interpret the deterministic realisation of the graphon, where each pair of spins interacts with a strength (vv ′ /L 2 ).Both interpretations lead to the same physics in the thermodynamic limit-this follows directly from theorem 2. From Eq. ( 7) we derive (see SM) σ z (x) = tanh(βcx) and M z = ln(cosh(c)) c where c is the real-valued solution of the equation and PL 2 is the PolyLogarithm function of order 2. The critical inverse temperature β c is β c = 3: the supremum of the LHS of the above equation for c ∈ [0, ∞].
In Fig. 1 we plot the total magnetisation M z and the local magnetisation σ z (x) versus β based on our analytical solution.We also perform finite-size Monte-Carlo numerics for M z for increasing system size (by constructing stochastic realisations of W (x, y)) and demonstrate convergence to our analytical solution.We compare these results to the graphon W (x, y) = 1 4 .As the temperature increases both systems undergo a second-order phase transition characterised by typical mean-field exponents.For the graphon W (x, y) = xy, however, the convergence to a fully ferromagnetic state at zero temperature is slower.This convergence can be determined analytically by expanding Eq. ( 8) for large β and substituting into β .There is thus a direct linear convergence of the magnetisation to unity with temperature T = 1 β versus the exponentially fast convergence associated with the homogoneous W (x, y) = const.case.
This slow convergence is a result of the 'left boundary' of the system.In Fig. 1b) we see that the local magnetisation at small values of x, where the spin-spin couplings are very weak, is very small even deep in the ferromagnetic regime.This 'boundary effect' means the T = 0 state has a finite magnetic susceptibility to changes in temperature, i.e. dM z dT | T =0 = −2 ln(2).In the homogeneous case we have dM z dT | T =0 = 0. Whilst both systems are mean-field in terms of their universal behaviour, they exhibit very different physics in the ferromagnetic regime.
Transverse field Ising model -We now consider a quantum example: the transverse field Ising model.Our integral equation is (setting J x = J y = h y = h z = 0 and h x = −h, J z = −1 in Eq. ( 1) dy.
(9) We focus on the ground-state by taking the limit β → ∞.We can again use Banach's fixed-point theorem here to prove the existence of a disordered-ordered phase transition with an upper bound of the critical field strength h c given by the supremum of the marginal of the graphon.
We now consider some specific examples.First, taking the Erdős-Rényi graphon W (x, y) = p straightforwardly yields the solution M z = 1 − h 2 p 2 consistent with a rescaled TFI model with all-to-all coupling [28].
There are, however, other, less trivial graphons for which an exact analytical solution for the ground state properties can be found.Consider the separable case W (x, y) = √ xy.Some algebra on Eq. ( 9) (see SM) leads Integrating the expression (see SM) for the transverse magnetisation then gives the following closed form for the total transverse magnetisation density The total longitudinal magnetisation density can also be obtained in closed form (see SM).Our methodology has yielded an analytic expression for the magnetisation (in the thermodynamic limit) of the transverse field Ising model on a complex, highly inhomogeneous graph structure.
In Fig. 2 we plot these solutions alongside those for the constant graphon.The left boundary of the system, which has very weak z − z coupling, modifies the physics of the system and makes it more susceptible to the transverse field than the all-to-all case.The transverse-field susceptibility versus site-index x can be derived from Eq. ( 10) yielding lim h→0 , where δ(x) is the Dirac-delta function.There is a singularity in the susceptibility on the left boundary of the system at 0 fieldstrength in the ferromagnetic regime.This is not present in the all-to-all model.Critical exponents for the magnetisations at the phase transition can be found via expansion of the analytical results and these are consistent with the mean-field universality class and equivalent for the two graphons.
In Fig. 2 we also provide finite-size simulations of the ground state on random-exchange realisations of W (x, y) using Density Matrix Renormalization Group (DMRG) [30] calculations on a Matrix Product State ansatz.We reach system sizes on the order of ∼ 100 spins, observing convergence to our analytical solution.We verify this convergence for local observables and non-local ones, where exact, analytical results can be obtained via the mean-field approximation.
Importantly, from these tensor network numerics we can go beyond mean-field theory and obtain the entanglement entropy of the ground state on a finite systemsomething currently inaccessible to our continuous formalism.This is non-zero and diverges logarithmically with the transverse field strength as criticality is approached: h → 0.5 − .We also find the entanglement only depends on the ratio x = N/L, where N is the partition size.This scaling is reminiscent of the entanglement properties of the all-to-all transverse field Ising model [23].Here we observe it in a heterogeneous dense graph system, suggesting a possible universal mechanism underpinning the scaling of entanglement entropy in these models.
Conclusion -We have successfully utilised tools from graph theory to derive a set of integral equations which describe the physics of generic spin models with a large density of interactions in the thermodynamic limitwhether classical or quantum.Our formalism straightforwardly reproduces known results and, most importantly, can be used to uncover the equilibrium properties of more complex systems.We observe how inhomogeneity in the underlying graphs alters the magnetic properties of the system.
Our work opens a up a number of further avenues for future research.Firstly, extending our formalism to describe the out-of-equilibrium dynamics of a spin system on a dense graph is a natural direction.Whilst an analytical solution is known for the all-to-all case (W (x, y) = 1) on the Lipkin-Meshkov-Glick model (a model whose dynamics was recently realised on a quantum simulator [24]), our graph-theoretic approach could open up solutions for a whole range of dense graphs.The quantum fluctuations which deviate finite-size results from the mean-field case would be stronger here.
Secondly, graphon estimation is the process of estimating the continuous graphon W (x, y) from which a given finite graph G could have been drawn from [31][32][33].Therefore when studying spin models on a large, connected structure (the structure need not necessarily be dense, graphon estimation can be done for quasi-sparse graphs too [32,33]) one can estimate the graphon W (x, y) and solve our equations to obtain an approximate solution to the equilibrium physics of the system.We restate the Hamiltonian from the main text ) where all definitions are retained and α = x, y, z.We now prove Theorems 1 and 2 from the main text, which are restated below.
..) be a sequence of finite-size graphs generated as stochastic or weighted realisations of the graphon W (x, y).Then for a given inverse temperature T = 1/β the macroscopic properties of the equilibrium states of the sequence of Hamiltonians H(G L ) L∈N = (H(G 1 ), H(G 2 ), ...) converges and are determined by the solution of the following coupled integral equations L and G W L be the stochastic and weighted realisations on L vertices of a graphon W (x, y) respectively.For an arbitrary set of real, finite values for the parameters {J x , J y , J z , h x , h y , h z } the following is true which vanishes in the limit L → ∞.
We will first prove Theorem 1 by assuming that Theorem 2 is true.Then we will prove Theorem 2 to complete the proof.
We first perform a mean-field treatment of H(G L ) for some arbitrary graph G L with adjacency matrix elements A v,v ′ and take L → ∞.Let σα v = ⟨σ α v ⟩ + δα v , substitute it into the Hamiltonian and ignore terms of order δ2 .The result is (up to a constant): (S4) Within this mean-field approximation the equilibrium state of the system is given by (S5) Where the reduced density matrix on each site ρ v is, explicitly (in the basis spanned by the eigenstates of σ z ), the following 2 × 2 matrix: where we have defined By taking the expectation values ⟨σ α v ⟩ associated with ρ v we find the λ α v variables must obey the following selfconsistency relation: with v = 1...L and α = x, y, z.The set of values {λ α v } with v = 1...L, α = x, y, z which solves the 3L non-linear equations described by Eq. (S7) thus fully characterise the mean-field equilibrium state associated with H. Now we wish to take the continuum limit of Eq. (S7).First, we define the following: We assume that the adjacency matrix has been generated as a weighted realisation of some graphon W (x, y), i.e.
Substituting this all into Eq.(S7) gives us Now we take L → ∞ which implies L 2 Iv×I v ′ W (x, y)dxdy → W (x, y) and the summation becomes an integral.We can then write down the coupled, continuous mean-field equations These govern our system in the thermodynamic limit of the sequence of graphs generated from the graphon W (x, y).Whilst we explcitly used the weighted realisation G W L of W (x, y), Theorem 2 tells us that the equilibrium properties of the system that arise as the solution of Eq. (S9) are equivalent for both G W L and G S L as L → ∞.Thus these equations govern the properties of any sequence of finite graphs which converge to W (x, y) -not just weighted ones.
The equations in Eq. (S9) are coupled, non-linear Fredholm integral equations with the graphon acting as the kernel.From the solution set {λ x (x), λ x (y), λ y (z)} to these equations we can obtain the on-site magnetisations via and the total magnetisation density is In order to complete the proof of Theorem 1 we need to prove Theorem 2 which was assumed at the end of the last section.We recall a Lemma proven in Ref. [28], which we will be helpful in completing the proof.We begin by defining the following operator: where A S ij are the matrix elements of G S L and A W ij are the matrix elements of G W L , finite stochastic and weighted realisations of some graphon W (x, y).
We proceed to evaluate the eigenvalues of the operator D α L .As such, consider its eigenstates |σ α 1 , ..., σ α L ⟩, where σα i |σ α 1 , ..., σ α L ⟩ = µ i |σ α 1 , ..., σ α L ⟩, with µ i = ±1, depending on whether the ith spin is pointing 'up' or 'down' in that basis.We will define µ ij := µ i µ j and consider the eigenvalue We we will proceed to show that, independent of the eigenstate |σ α 1 , ..., σ α L ⟩, this eigenvalue grows, at most, as L 1/2 in the large L limit.From there we can invoke Weyl's inequality to show that the eigenvalues of D L grow asymptotically as L 1/2 and subsequently invoke Lemma 1 (from above) to complete the proof of Theorem 1.
To prove the L 1/2 growth of Eq. (S12), we begin with the Hoeffding inequality.This states that for independent random variables Y 1 , ..., Y n for which ) with the factor of 2 stemming from the fact we have incorporated both the upper and lower Hoeffding bounds together.We will apply this bound to Eq. (S12).
Let us construct the set X := {X ij } which consists of the L(L−1) ) denotes the expected value.In our case, however, all of the X ij quantities are not independent-the sign of X ij is determined from the sign of X ik and X kj .This can be dealt with by applying the Hoeffding bound to the X ij 's with positive and negative sign separately.
A given eigenvector |σ α 1 , ..., σ α L ⟩ will consist of a number of spins pointing up (σ i = +1) and and the remainder pointing down (σ i = −1).Define M = L k=1 σ k , then it can be checked that of the L(L−1) are negative.We therefore partition the set X := {X ij } into two sets as follows We would like also to keep track of the values A W ij µ ij , so we partition the set We can now invoke the Hoeffding bound on each set separately, since each set now consists on independent random variables.We have two sums on which to invoke the bound: x A l and S B := x B l , where we use x A l and x B l to refer to elements from x A and x B respectively.Likewise we will use xA We combine these bounds to obtain: Since we fixed the magnetisation M of the eigenstate in deriving the above bound, we should take the union bound over eigenstates with magnetisation M .We will denote an eigenstate with magnetisation M as |σ M ⟩, resulting in the new bound Now observe that the following is true where This leads us to the following where the union bound has again been used.Now observe that if Lπ 2 L for large L, and also that using 2 L = e Lln (2) , we then arrive at which vanishes unless γ ≤ 3 2 .This leads us to Now from Weyl's inequality we know that the eigenvalues of the operator ).From here we can invoke Lemma 1 with κ = 1/2 and Theorem 2 is proven.

B. Appendix B: Analytical solution of the classical
Ising model for W (x, y) = xy We wish to solve Eq. ( 7) with W (x, y) = xy, i.e. identify the function λ z (x) which solves We can perform the integration here analytically.First, perform integration by parts and expand into exponentials: Now the integral on the RHS can be dealt with by observing that e −2y (1+e (S23) We can evaluate the series by splitting up the numerator and using known results, We can then use this result to reduce Eq.(S20) to Eq. ( 8) from the main text: We wish to solve Eq. ( 9) in the main text with W (x, y) = √ xy and β = ∞, i.e: again observing that this implies λ z (x) = √ xg(h) where g(h) is some real-valued function of h gives which we can solve for g(h) (we restrict h and g(h) to be positive real without loss of generality) by a series of substitutions.Yielding (S28) which has the solution as in the main text.The transverse and longitudinal magnetisations are determined by the integrals dx respectively.The first, by direct integration, yields Eq. ( 11) in the main text.The second can be done by an extensive series of trigonometric substitutions and results in the following closed form expression where s = (−1 + 3h) √ 1 + 6h.

D. Appendix C: Further example graphons
In this section we consider further graphons which were not treated in the main text but frequently appear in the literature on graphons.
Stochastic block models -The Stochastic Block Graphon is typically utilised in statistical analysis of networks because they are useful in uncovering clustering in networks [35].The graphon can be expressed as with p ij = p ji and the X i specifying disjoint sub-domains of [0, 1] such that ∪ k i=1 X i = [0, 1].We write ∆X i to indicate the width of the interval X i .The continuous mean-field equations then take on the following form Observe that we can immediately infer from this that λ α (x) is constant across each of the domains X i .We can thus define λ α i = λ α (x) ∀x ∈ X i and reduce Eq.(S33) to a series of equations which become increasingly complicated to solve as the number of clusters does.In the case of a single cluster we recover the case of an Erdős-Rényi graph.
Growing uniform attachment -The growing uniform attachment graphon is given by W (x, y) = 1 − max(x, y) [36].The graphs which are finite realisations of this graphon will consist of nodes in which the average connectivity of a node varies uniformly across the graph.Such graphs are therefore highly inhomogenous in their average vertex connectivity.Substituting W (x, y) = 1 − max(x, y) into Eq.(S9) and differentiating the left and right hand sides twice with respect to x leads us to the following coupled second-order ODEs (S34) with α = x, y, z, boundary conditions λ α (1) = 0 and Maximally irregular graph -The maximally irregular graph is the finite connected graph where each site (other than one pair) has a different co-ordination number to any others [28].Taking the thermodynamic limit of the adjacency matrix results in the graphon and the integral equations in Eq. (S9) reduce (upon differentiation) to the following three coupled first order ODEs with boundary conditions λ α (1) = 0 ∀α.Such equations are known as functional differential equations and have been studied extensively in both Mathematics and the applied Sciences [37].

E. Appendix D: Numerical details
Classical Ising model -For the finite-size data plotted in Figure 1 of the main text we used Monte-Carlo simulations.Specifically, for a given L we drew a finite randomexchange realisation of the graphon W (x, y) = xy and for a given temperature β utilised the Metropolis-Hastings algorithm to generate N Samples = 5000 for the Magnetisation Density M z .We used a Markov chain length of 250 between each sample and threw away the first 1000 samples.For each L we took 100 stochastic realisations of the graphon W (x, y) and averaged our results over this.There are thus two sources of statistical error in our simulations: the error from sampling a finite number of stochastic realisations and the error from taking a finite number of Monte-Carlo samples.In Fig. 1 we plot the standard error on the mean from both of these sources, the values are negligible in comparison to the scale (0 → 1) of Fig. 1 in the main text.
Transverse Ising model -For the data plotted in Figure 2 of the main text we used the Density-Matrix-Renormalisation-Group (DMRG) algorithm to find the ground-state of the transverse field Ising model.For a given L we drew a finite stochastic realisation of the graphon W (x, y) = √ xy.Then, for a given field strength h we took a random Matrix Product State with a small bond dimension χ and successively performed DMRG sweeps, letting the bond-dimension double every 4th sweep until the energy converges to within 0.1% of that for the previous bond dimension.There is thus only one source of statistical error in this simulation: the error from sampling a finite number (100) of stochastic realisations.In Fig. 1 we plot this error as a percentage and observe that it is on the order of 0.1%.The ordering of the sites (from left to right) of the Matrix Product State was taken to be identical to the ordering v = 1...L of the sites of the graph.

F. Appendix E: The Graphon as the Limit Object of Dense Graph Sequences
We provide mathematical details on how the graphon W is the limit object of a sequence of dense graphs (G n ) n∈N where n is the number of vertices.This Appendix closely follows Ref. [38], although the theory on graph limits was first developed in Ref. [17].The interested reader should consult either of these for more detail.
Consider two simple graphs F and G, where we define the number of vertices of F to be k and that of G to be n.A homomorphism from F to G is a map which preserves edges.This means that given an edge (i, j) ∈ E(F ) (here E(F ) is the edge set of F ), and a homomorphism h, there is always an edge (h(i), h(j)) ∈ E(G)-the set of edges of G. Let hom(F, G) indicate the number of homomorphisms from F into G.The homomorphism density t(F, G) is then defined to be The homomorphism density is the probability of a random map from the graph F to the graph G being a homomorphism, since n k is the total number of maps from a graph with k vertices to a graph with n vertices.Suppose that instead we are given a graphon, such as W G -the stepped graphon corresponding to the graph G which is defined as n] (with A being the adjacency matrix of G).In this case, the homomorphism density is Supplementary Figure 1.a) Standard error on the mean for the Monte-Carlo calculations of M z for the classical Ising Model on stochastic finite-size realisations of W (x, y) = xy.Top plot is standard error on the mean from 5000 Monte-Carlo samples of M z at a given β and L. Data points are averaged over 100 stochastic realisations of W (x, y).Bottom plot is the standard error on the mean from the 100 stochastic realisations of W (x, y), averaged over the 5000 samples taken for each realisation.b) Relative standard error on the mean (standard error on the mean as a percentage of the mean) for the ground state energy of the transverse field Ising model calculated via DMRG.Standard error is that originating from the 100 stochastic realisations of W (x, y) = √ xy at a given h and system size L.
defined to be Here the same definition holds for any arbitrary graphon W .
The homomorphism density with reference to a finite graph F indicates the relative likelihood of the graph G or more generally graphon W containing an instance of F inside of it.If two graphs or graphons have similar homomorphism densities for all simple graphs F , then these graphs are similar.The definition of convergence of a sequence of graphs hinges precisely on this concept.Definition 1 (Convergent Graph Sequence) A sequence (G n ) of simple graphs with V (G n ) → ∞ as n → ∞ converges if the subgraph densities t(F, G n ) converge for all simple graphs F .
The above definition gives allows us to precisely define in what sense W can be considered a limit object.
Theorem 3 (Lovasz, 2012 [38]) Let (G n ) be a sequence of simple graphs with V (G n ) → ∞.If (G n ) converges, there exists a graphon W such that t(F, G n ) → t(F, W ) for all simple graphs F .
The above theorem tells us that if the sequence (G n ) converges, then there exists some limit object-the graphon-which captures the limiting homomorphism density counts of the sequence of graphs for all simple graphs.
There is a second, equivalent, definition of convergence which us allows us to define W as an appropriate limit of a sequence of dense graphs.This definition utilises the cut distance of two graphs.where the infimum is taken over all vertex re-labelings ϕ of W and ψ of W ′ .The supremum is taken over all measurable subsets S and T of [0, 1].
The cut distance is a metric on the space of graphons (up to weak isomorphism).It maximises the difference between the integral of the two graphons on measurable intervals S and T which together form a box S × T .This step can be thought of as maximising the difference in edges between those vertices contained in S × T .The infimum is then taken on that chosen interval over all measure preserving maps, in order to ensure that the cut distance is zero for weakly isomorphic graphons.The following theorem can then be proven from the above definitions.
This theorem provides alternative definition for the graphon as a limit object.In this definition, we envisage instead the pixelated adjacency matrix of the sequence of simple graphs G n approaching (via the cut distance) that of the limit object W .
Importantly, the above definitions and theorems can be generalised to sequences of weighted graphs by requiring the graphs to have uniformly bounded edgeweights.Moreover, we emphasise that these limits only make sense for sequences of dense graphs, because it can be shown that sparse graph sequences always have as their limit the graphon W (x, y) = 0 for all x and y.

Theorem 2
Let f (H) = − 1 Lβ ln(Tr(exp(−βH))) be the free energy density of a d L × d L many-body Hamiltonian, with β ∈ R ≥0 and d the dimension of the local Hilbert space.Let G S

FIG. 1 .
FIG. 1. Magnetisation of the classical Ising model for the graphons: W (x, y) = xy and W (x, y) = 1 4 .a) Total Magnetisation density M z versus inverse temperature β for L → ∞.Black dashed-dotted lines give the asymptotic derived by taking the large β limit of the respective closed-form equations.Red circles correspond to Monte Carlo simulations of finite, L = 800, randomly sampled graphs G S L derived from W (x, y) = xy.Bottom) Percentage difference in M z for the exact result in the thermodynamic limit versus finite-size Monte Carlo simulations at several L (crosses are L = 100, triangles L = 200, squares L = 400 and circles L = 800).For each β and L, 100 stochastic samples G S L are realised and the data (both top and bottom plots) is averaged over these.Further details are provided in the SM.b) On-site magnetisation σ z (x) versus β and x for the graphon W (x, y) = xy in the thermodynamic limit.

FIG. 2 .
FIG. 2. Properties of the ground state of the transverse field Ising model on the graphon W (x, y) = √ xy.Results for the constant graphon are included for reference.a) Energy density versus transverse field strength h.Orange line represents the analytic solution in the thermodynamic limit.Markers represent numerical calculations averaged over 100 finite stochastic realisations of W (x, y) = √ xy on L = 400 sites.Inset) Percentage difference between the ground state energy calculated on L = 100, 200 and 400 (cross, triangle and circle marker respectively) site random-exchange realisations of W (x, y) and the exact solution for L → ∞. b) Total transverse (unfaded) and longitudinal (faded) Magnetisation densities of the ground state.Inset) Two-point correlator ⟨σ x ( 1 4 )σ x ( 3 4 )⟩.c) Von-Neumann Entanglement Entropy (EE) of W (x, y) = √ xy averaged over 100 stochastic realisations on L = 400 sites.The partition is between the first xL sites of the system and the remaining (1 − x)L sites.The red curve corresponds to the entanglement entropy at x = 2/3 for h = 0 → 1. Dotted black line is the fit EE(x = 2/3) = −0.136log 2 (h−0.5)−0.089.d) Analytical result for the on-site magnetisation σ x (x) versus transverse field strength h and position x for the graphon W (x, y) = √ xy.
A: Proof of Theorems 1 and 2.

Theorem 2
2 and the three functions λ α (x), with α ∈ {x, y, z}, each being continuous, real-valued and defined over the domain[0, 1].Let f (H) = − 1Lβ ln(Tr(exp(−βH))) be the free energy density of a d L × d L many-body Hamiltonian, with β ∈ R ≥0 and d the dimension of the local Hilbert space.Let G S

Lemma 1
Let (A L ) L∈N = (A 1 , A 2 , ...) and (B L ) L∈N = (B 1 , B 2 , ...) be two sequences of many-body Hermitian matrices.The matrices A L and B L in the sequence have size d L × d L , with d fixed and the dimension of the local Hilbert space.Let D L = A L − B L and λ D Max be the largest (in terms of the absolute value) eigenvalue of D L .If |λ D M ax | = O(L κ ) then |f (A L ) − f (B L )| = O(L κ−1 ), which vanishes for κ < 1 as L → ∞.

4
are positive and L 2 −M 2 4

M 2
i and xB i to refer to individual elements of xA and xB respectively.Then observe that, due to E(X ij ) = A W ij µ ij we have E(S A ) = ) then gives two bounds:

1 0
tanh(cx)dx = ln(cosh(c)) c .C. Analytical solution of the ground-state of the transverse field Ising model for W (x, y) = √ xy