key: cord-0893709-k3hfd8j9
authors: 'Odor, G'eza
title: Non-universal power-law dynamics of SIR models on hierarchical modular networks
date: 2021-03-12
journal: Physical review. E
DOI: 10.1103/physreve.103.062112
sha: 9b4b2a3dcb567e19601334e9287bb588ef6033d9
doc_id: 893709
cord_uid: k3hfd8j9

Power-law (PL) time dependent infection growth has been reported in many COVID-19 statistics. In simple SIR models the number of infections grows at the outbreak as $I(t) propto t^{d-1}$ on $d$-dimensional Euclidean lattices in the endemic phase or follow a slower universal PL at the critical point, until finite sizes cause immunity and a crossover to an exponential decay. Heterogeneity may alter the dynamics of spreading models, spatially inhomogeneous infection rates can cause slower decays, posing a threat of a long recovery from a pandemic. COVID-19 statistics have also provided epidemic size distributions with PL tails in several countries. Here I investigate SIR like models on hierarchical modular networks, embedded in 2d lattices with the addition of long-range links. I show that if the topological dimension of the network is finite, average degree dependent PL growth of prevalence emerges. Supercritically the same exponents as of regular graphs occurs, but the topological disorder alters the critical behavior. This is also true for the epidemic size distributions. Mobility of individuals does not affect the form of the scaling behavior, except for the $d=2$ lattice, but increases the magnitude of the epidemic. The addition of a super-spreader hot-spot also does not change the growth exponent and the exponential decay in the herd immunity regime.

Human infectious diseases usually start with an exponential growth, as the number of healthy neighbors is high, thanks to the small world connectedness of societies. Thus a full, infinite dimensional graph approximation, described by mean-field behavior is valid. In finite d dimensions this evolution is slower and in case of the Susceptible Infected Recovered (SIR) process [1] , the simplest model for epidemics with immunization, the actual number of infected individuals follows a scaling behavior to leading order: I(t) ∝ t d−1 . However, this is true in the supercritical phase, where the reproduction number of epidemiology is R 0 > 1. By reducing R 0 a continuous phase transition to a non-endemic state happens and right at the critical point: R 0 = 1 the number of infected individuals grows algebraically, characterized by the so-called initial slip exponent η in the scaling law I(t) ∝ t η of statistical physics [2] [3] [4] . For different models the value of η is known in different Euclidean dimensions of the substrate graphs [4, 5] , in homogeneous systems.

In the case of quasi-static, or quenched heterogeneity much less is known. According to the Harris criterion [6] , the disorder is irrelevant for the critical behavior of SIR, belonging to the Dynamical Isotropic Percolation (DIP) universality class [4, 7] . This means that the correlation length exponent ν ⊥ fulfills the inequality: dν ⊥ > 2 and weak disorder decreases under coarse gaining becoming unimportant on large length scales. But this criterion is a necessary and not sufficient condition for the stability of the impure fixed point and has been found to fail in certain models of nonequilibrium statistical physics [8] .

Containment measures can push R 0 below 1 by lowering the infection rate or the graph dimensions for example. Large amount of Covid data have provided various growths functions of I(t) with PL, exponential or mixed time dependence [9] in different countries and at different times.

Very recently a basic SIR like model has been investigated on Euclidean lattices with inhomogeneous infection rates, including the possibility of mobility [10] . A striking numerical conclusion was drawn, that in the presence of super-spreader hot spots, where the infection rate is much higher than the average, the epidemic does not vanish exponentially fast by herd immunity, but decays in a slow PL manner way. This behavior was paralleled with the Griffihts Phase (GP) phenomena [8] , which occurs near the critical points of phase transitions in strongly heterogeneous systems. At first glance this seems surprising, because for having GP one wold need long surviving rare regions (RR), in which the activity disappears exponentially slowly by the region size: τ t ∼ e V , while in the SIR model recovered individuals cannot be re-activated but become inactive forever. Thus the effective topological dimension of an infected region decreases quickly as herd immunity develops. The authors of Ref. [10] suggest a clue for this strange behavior by the interplay of SIR processes and diffusion. Indeed, mobility can increase the effective dimension of systems and in the infinite diffusion limit mean-field behavior emerges, with a PL decay at criticality. Furthermore, mobility of individuals can transform recovered sites susceptible again and Susceptible Infected Susceptible (SIS) type of scaling behavior, belonging to the Directed Percolation (DP) [2] [3] [4] universality class may be observed.

The question is whether the diffusion is strong enough to counterbalance the spontaneous reduction of the effective dimension caused by the recoveries in a finite system. In Ref. [10] no systematic investigation has been provided to understand this better. Furthermore, modern human societies cannot be described by regular Euclidean lat-tices, rather by small world graphs, restricted by containment measures in the course of defense. Therefore the infinitely strong diffusion limit is unrealistic. To describe epidemics usually meta-population models are used, built from internally strongly connected modules, which are interconnected via diffusion of scale-free graphs [11] .

Here I advance another assumption for modeling societies with lock downs. This is to be done using hierarchical modular networks (HMN), where the modules can be families, villages, towns, countries and continents with random intra-module connections, while the modules, embedded in the 2d space are interconnected via long links with geometrical distance decreasing probabilities. In this work I investigate the dynamical behavior of SIR like models on such HMN-s, considering hot-spot heterogeneity as well.

Another very recent publication [12] claims that empirically the average fraction of infected people decays over time algebraically and tries to understand it via SIR like models on fixed graphs. The authors conclude that even non-Markovian description fails to explain the empirical data and conjecture that time-varying, human contact graphs are needed to produce slow dynamics. Indeed, such graphs may be more realistic, even if containment measures freeze the mobility. In this work I also present results for the mobility effects on the fixed graph SIR model results.

In this section I describe the HMN networks used for the simulations. The network generation starts at the top level, by connecting neighbors to the N nodes via edges with the probability p 0 = b( 1 2 s ) lmax . Here s determines the type of the network and b is a control parameter. Then, further random long links are added by level-to-level from top to bottom, similarly as in [13, 14] , excluding self-connections. The levels: l = 0, 1, ..., l max are numbered from bottom to top. The size of domains, i.e. the number of nodes in a level, grows as N l = 4 l+1 in case of the 4-module construction, related to tiling of the 2d base lattice (see inset of Fig. 1 ). The probability of random, intra-module links at level l is

With this construction he average degree k of nodes is related to b as shown in Table I . Nodes are connected in a hierarchical modular way as if they were embedded in a regular, two-dimensional lattice (HMN2d) as shown by the adjacency matrix on Fig. 1 , similarly as in [13, 14] . The 4 nodes of the level l = 0 are set to be fully connected. These HMN-s, possess increasing edge density from top to bottom levels. Such topology has been shown to be suitable to describe activity localization and for the emergence of potential RR effects [13, 15] . One can make a correspondence with spatially embedded networks of type discussed in [16] . These networks have long links, with algebraically decaying probabilities in the Euclidean distance R as

Single connectedness of networks is not required, a typical l max = 5, s = 4 sample with N = 4096 nodes and 37240 edges contains 3 strongly and 1 weakly connected components. Other, randomly selected networks showed statistics and invariants within a few percent difference. The modularity coefficient of the networks is high: Q > 0.94, defined by

where A ij is the adjacency matrix and δ(i, j) is the Kronecker delta function. The Watts-Strogatz clustering coefficient [17] of a network of N nodes is

where n i denotes the number of direct edges interconnecting the k i nearest neighbors of node i. In randomly generated HMN2d-s with N = 4096 this is roughly C = 0.31, which is more than 122 times higher than that of a random network of same size C r = 0.002548, defined by C r = k /N . The average shortest path length is defined

where d(i, j) is the graph distance between vertices i and j. For several typical networks L = 10.44 is found, which is larger than that of the random network of same size: L r = 4, computed from the formula [18] :

Hence these are small-world networks, according to the definition of the coefficient [19] :

because σ ≃ 47 is much larger than unity. We can estimate the effective topological (graph) dimension d T , using the breadth-first search (BFS) algorithm, defined as

by counting the number of nodes N (r ′ ) with chemical distance r ′ ≤ r within a large sample averege of trials started from randomly selected seeds. The dimension d is estimated for different l max , b and s values. For s = 3 we can't find the true d T → ∞ due to the finite size cutoff. As one can see on The average number of nearest neighbors measured by the BFS algorithm as the function of graph distances from randomly selected initial seeds. Different b and s values and sizes lmax were investigated as shown by the legends. Lines correspond to least squared PL fits for 10 < r < rcut, where rcut was estimated visually.

Time dependent simulations were performed from randomly selected infected seed initial conditions. This means that at t = 0 randomly selected i 0 number of nodes were set to the infected state: x(i) = 1 in an otherwise fully susceptible system x(j) = 0 for j ∈ (1, N ). Usually i 0 = 1 was used, to describe epidemics from single sources, but multiple source cases have also been considered. This source triggers epidemic avalanches, used in statistical physics to investigate the so called critical initial slip phenomena [3] . The HMN2d graphs considered for extended simulations has l max = 6, 7, 8 levels, containing N = 16384, 65536, 262144 nodes, respectively. For testing purposes d = 2 and d = 3 dimensional Euclidean lattices with linear sizes L = 100, 1000, 2000 and periodic boundary conditions were also investigated. At times t = 1, 2, 3, ..., t max vector elements of the updated state variables are set to x ′ (j) = 1, with probability λ, provided they were in susceptible state before x(j) = 0 and had any infected neighbors. Infected nodes recover to the state x ′ (j) = −1 with probability ν. To study mobility effects the code also performs an exchange of states with that of a randomly selected neighboring node, if the reaction conditions described above are not satisfied. Following a full sweep of nodes the old state vector is updated with the new one : x(j) = x ′ (j), corresponding to one Monte Carlo step (MCs). Throughout this study I measure time in MCs units. Thus, Stochastic Cellular Automaton (SCA) like updates have been used without the loss of generality.

To prove equivalence of the scaling behavior with that of the SIR model the simulations were tested on d = 2, 3 dimensional lattices, with periodic boundary conditions, for which the DIP universal scaling exponents are tabulated [4] . The density of infected nodes I m (t) = 1/N N i=1 δ(x i , 1) is measured in each sample run m and the spatio-temporal size S m = N i=1 T t=1 δ(x i , 1) of the avalanches is calculated, where T denotes the maximal duration of the epidemic avalanche. Averaging over I m (t) of the independent samples we get I(t) = I m (t) . To determine PDF of S m a histogramming algorithm is used on the results of thousands of realizations, started from random initial conditions. That means random initial infected site locations, as well as random initial graph configurations in the case of the HMN2d networks. At the critical point the PL behavior of the PDF decay tail defines the exponent τ as : p(S) ∝ S −τ . Furthermore, the avalache survival probability P (t) is also determined, which scales at the critical point as: P ∝ t −δ . To obtain more precise exponent estimates the local slopes of I(t) is deretmined by

using t − t ′ = 4. Similarly one can analyze the P (t) results. These effective exponent curves, plotted as the function of 1/t veer up or veer down super-or subcritically. At the critical point no curvature is expected in case of simple PL sub-leading corrections to scaling and on can read off η by extrapolating to 1/t → 0 on the vertical axis. In d = 2 dimensions square lattices of linear size: L = 1000, 2000 were used and the critical scaling results were found to be in full agreement with those of the Dynamical Isotropic Percolation (DIP) class [4, 20] . At the critical point λ c = 0.4059(1), determined by the local slope analysis, I estimate: η = 0.59(1) as compared to η DIP 2d = 0.586 [20] (see Fig.3 ). In d = 3 cubes of linear sizes: L = 100, 160 were used and at the critical point λ c = 0.2198(2) the growth exponent η = 0.53 (2) was found in comparison with literature value η DIP 3d = 0.536 [20] . Above the critical point we can observe the expected scaling, characterized by the exponent η = d − 1, before the size cutoff turns on. However, we can also see nontrivial corrections, causing an overshoot of the effective exponents if λ is slightly above λ c . This correction is more pronounced in the d = 3 case, where smaller sized lattices could be accessed than in d = 2 dimensions.

In case of the HMN2d graphs as the first step I determined the growth behavior of I(t) at λ = ν = 1 for different s values. As one can observe on Fig. 4 , PL-s seem to occur for s = 4, where the topological dimension of the graphs is finite. For s = 3 the dimension is infinite and we can see faster then PL growth behavior. For s > 4 the spatial dimension vanishes: lim N →∞ d T → 0 and we can find slower than algebraic initial growth of I(t). Note, that for obtaining collapse of densities with different sizes the l max = 8 data was multiplied by a factor of four. I have also investigated the initial i 0 de- pendence by varying it from i 0 = 1 to i 0 = 16, because a recent study of the mean-field SIR model [21] suggested the possibility of i 0 dependent critical exponents. The present simulations show that the number of initial seeds scales up the magnitude and the duration of the I(t) curve, but the slope of the PL does not change.

Next I concentrated on the s = 4, b = 1 case and determined the infection probability dependence as shown on Fig. 5 using ν = 1 − λ. In this case one initially thinks of λ dependent initial PL-s by looking at the results. When we calculate the local slopes of I(t) we can observe a phase transition point at λ c ≃ 0.3 characterized by η = 1.4(1) as it appears in the inset of Fig. 5 . This exponent is much bigger than that of the DIP universality class value for the d = 3 case: η DIP 3d = 0.536(10) [5] . Since η DIP decreases further, when the spatial dimension increases the HMN2d exponent is in conflict with the expectation of a homogeneous system projection for d T = 3.5, obtained by BFS analysis for this network. Thus, the topological disorder alters the critical point scaling behavior.

We can see that the supercritical curves with λ c < λ < 0.45 veer up, while those with λ > 0.45 veer down, converging to η s ≃ 2.5(10), in agreement with the result obtained by the BFS dimension measurements, according to which we expect η s = d T − 1 ≃ 2.5. Note, that the small variation of the effective exponents in the narrow scaling region before the finite size cutoff may suggest the wrong conclusion of λ dependent scaling exponents. We have no reason to believe in such non-universal exponents here in the lack of long surviving RR-s. This will be more obvious later, at the b = 0.4 case, where the exponents of the two fixed point, the critical and the supercritical one are more distant. I have also performed test runs by starting the system from homogeneous, 80% randomly infected initial states near criticality and found simple exponential decays of I(t), ruling out of a GP like behavior. Note that the small zig-zags at t ≃ 20 are numerical artifacts, coming from the combination of SCA updating on HMN2d lattices. They do not appear in the regular lattice simulations. I have repeated these simulations for other b values at s = 4 and found similar results. For example at b = 0.4 the average degree is k = 6.3 and the topological dimension is d T ≃ 3. The curves of the effective exponents on Fig. 6 veer down for λ < 0.48 and veer up for λ > 0.48 as t → ∞, but finite size terminates the epidemic and causes an exponential cutoff. At λ c ≃ 0.475(1) we can read off an asymptotic PL growth behavior characterized by the exponent η = 0.8(1), much larger than the homogeneous, 3 dimensional system value: η DIP 3d = 0.53(2). One can also see clearly, that in the supercritical phase the effective exponent curves do not level off, but tend to η ≃ 2 = d T − 1, following the overshoot correction region. For other b values we can obtain similar results, the critical and the supercritical η exponents change continuously as shown in Table I , summarizing the results.

The total epidemic size statistics, determined for different countries, also show PL distributions [22] and a snowball model on a two-level, heterogeneous system was (2d, 3d) , or by the value of b for HMN2d. "+D" denotes the diffusive case, "+H" means the application of a single super-spreader hot-spot. suggested to describe it. Another, earlier, brain motivated study of a SIR like system, applied on HMN2d graphs also concluded PL-s for the fractional component sizes [23] using normalization and simulation methods. I have determined the PDF-s of the infection spatiotemporal avalanche sizes of the s = 4 case for different b and λ values near criticality. These distributions seem to exhibit PL tails before a bump at the end in case of supercriticality, corresponding to a giant component. Again, first I calculated these distributions for testing in case of Euclidean lattices and found results in good agreement with the corresponding critical DIP classes. In particular, in d = 2 I obtained: τ = 1.06(1) as compared to τ DIP 2d = 96/91 ≃ 1.055 of the DIP class [20] . In d = 3 this model provides τ = 1.20(2) with respect to τ DIP 3d = 1.188 of the DIP class [20] . For the supercritical λ-s in infinite systems the epidemic never stops, thus the p(s) distribution is singular: p(S) ∝ S −1 . Fig. 7 shows the results for HMN2d-s with s = 4 and b = 0.4, corresponding to average degree k = 6.3 and graph dimension d T ≃ 3. A PL fit for the intermediate "tail" region, which is s > 10 before the bump, corresponding to a giant component, results in: τ = 1.05(5), a slightly smaller exponent than that of the regular lattice: τ DIP 3d = 1.20 (2) . Several bumps can be seen, due to modules of different sizes. These log-periodic oscillations, superimposed on the PL-s are the consequence of the discrete scale invariance of the HMN2d graph and increase the numerical uncertainty of the fitting procedure. 

To simulate time-varying, human contact graphs I repeated the aforementioned spreading analysis for 2d lattices as well as for the HMN2d-s with s = 4, by allowing the diffusion of states. This means that SIR individuals are not fixed and can have different neighbors. The simulation program emulate this by an additional state exchange of x(i) with the state of a randomly selected neighbor x(j), provided the reaction conditions are not satisfied.

For the 2d lattice the supercritical scaling is invariant, but the critical point increases slightly to λ c = 0.413(1) and the growth exponents decreases to η = 0.25(2) (see Fig.8 ), deviating considerably from the 2d SIR exponent. This is the consequence of the site reinfecibility, the long term memory of sites is lost, causing SIS type of critical behavior. This exponent is close to the 2d DP universality class value: η 2d,DP = 0.2295(10) [2] [3] [4] . It is hard to determine it very precisely, due to the fi-nite size cutoff, but increasing the size from L = 1000 to L = 2000 just above λ c one can observe a up-bend curvature before the cutoff, which moves the estimate towards η 2d,DP . However, for the survival probabilty exponent we can obtain an estimate: δ = 0.20(1) before the exponential break down, which is far away from that of the DP value δ 2d,DP = 0.4505 (10) and also from the DIP value δ 2d,DIP = 0.092 [2] [3] [4] . Interestingly the diffusion increases the inactive phase, susceptible individuals can diffuse back, behind the growing epidemic front, becoming target of re-infection. The size distribution exponent τ does not change, neither the exponential decay in the herd immunity phase. In case of the HMN2d networks the diffusion does not seems to change the dynamical exponents. As we can see on Fig.9 , in the supercritical phase at λ = 0.5, the initial scaling behavior remains the same, characterized by η ≃ d T − 1 = 2.5, but the epidemic grows further, achieving a larger maximum value than in the frozen case. The inset of The inset of Fig. 9 shows the effective growth exponent results for different λ-s at b = 1 near the critical point. Again a critical point at λ c = 0.240(3) appears, lower than that of the frozen case (0.310(5)) and in supercritical phase η ef f estimates tend to η ≃ 2.5.

The scaling behavior of the epidemic size distributions P (S) exhibit insenitivity to the mobility as shown by main plot of Fig. 10 ), but the oscillations are even more pronounced. The fitted exponent for the tails before the bump, corresponding to the giant component τ = 1. 10(8) is in good agreement with that of the frozen network case.

In contrast with the 2d lattice, the mobility does not seem to alter the critical point scaling, probably because the strongly connected network structure does not allow recovery, keeping the long-time local memory of sites. C. Super-spreader hot spot in the presence of mobility I have also investigated the effect of a super-spreader hot spot in the presence of mobility as suggested in Ref. [10] . This was achieved by increasing the infection probability at a single site i = 100 to λ = 1. However, this change alone did not cause measurable difference in the epidemic, neither in the initial regime, nor in the size distributions in case of HMN2d-s. Fig. 11 shows the I(t) for s = 4 and b = 1 situation. In the herd immunity regime the decay remains exponential. Thus in the HMN2d graphs, where the epidemic propagation is restricted by the "containment" of modules neither the mobility nor hot-spots can lead to slow decay dynamics as in Ref. [10] . In case of 2d diffusive lattices at λ = 0.5 this single hotspot increased the size of I(t), by growing the maximum by 20%, as one can observe in the inset of Fig. 10 , but the decay remains exponential. Without diffusion, quenched hot-spots do not have strong effects in 2d, even a 10% concentration of randomly distributed ones cause similar growth with respect to the homogenous lattice and the peak time seems to decrease. The details of quenched disorder runs are presented in the Appendix.

SCA models following SIR rules on regular lattices and on hierarchical modular graphs have been investigated by numerical simulations. Scaling behavior at the epidemic outbreak was found in systems with finite topo-logical dimensions. The dynamical percolation behavior on d = 2, 3 dimensional lattices was fully confirmed and a nontrivial correction to scaling in the supercritical phase has been explored in detail.

To model human society with containment measures hierarchical modular graphs, embedded in 2d space (HMN2d) were used. A special set of graphs were considered, in which long-range connections decay with the geometrical distance in a PL manner and the topological dimension, together with the average degree can be tuned via a single parameter. I found that the critical behavior is altered by the topological heterogeneity, such that k dependent exponents occur for the number of newly infected agents as well as for the total spatio-temporal sizes of pandemics, which can be regarded avalanches, triggered by a single infected site. By changing continuously k , non-universal PL behavior emerges for the growth and epidemic sizes, characterized by the continuously varying exponents η and τ . Comparing exponents of a d T ≃ 3 HMN2d with those of a d = 3 Euclidean lattice provides different exponents, thus it turns out that the topological heterogeneity is relevant for modifying the dynamical scaling behavior.

On the other hand, the supercritical scaling is insensible to the heterogeneity, one can find the same growth laws as in the corresponding d dimensional Euclidean lattices. The epidemic size distributions also exhibit PL tail region at criticality, with exponents smaller, but close to those of the homogeneous lattices. They increase slightly with k .

Note, that as the scaling region, before herd immunity sets in, is narrow in finite systems and one can misleadingly think of λ dependent exponents for a given topology. Empirical data may also suggest this. Only a detailed, local slope scaling analysis, showing the scaling corrections allows one to determine the true asymptotic behavior. Smaller system sizes make the outbreak scaling region narrower as well as the maximum of I(t) smaller, but do not change the scaling exponents. Multiple sources also do not affect the exponents, but the scale up sizes of the epidemics.

Comparing values of η and τ with those from COVID-19 statistics we can see that a proper HMN2d model description would require larger k and graph dimensions to obtain agreement with real data. Fitting for the confirmed COVID-19 case statistics τ = 1.14−1.5 exponents are reported in [22] . The largest critical value investigated here was τ = 1.12(5) but this can grow further by increasing the control parameter b. Of course the epidemics are not necessarily critical, nor their connection network is finite dimensional. Furthermore, in reality the parameters change in time, thus for example the catastrophic case with a giant component can be avoided.

Effects of mobility and the possibility of hot-spots with large local infection rates have also been investigated. Diffusion washes out long term memory of sites and in the 2d lattice and we can observe DP universality class like initial slip scaling exponent η as in case of the memoryless SIS model. However, the avalanche survival exponent δ is different. Further investigation to clarify if this is a new universality class or diffusion strength dependence occurs would be needed.

In case of HMN2d the mobility did not change the late time exponential decay dynamics. It increases the size of epidemics, but does not seem to alter the critical exponents. Therefore, the HMN2d modular topology is efficient to keep the functional form of the frozen SIR epidemic dynamics. The addition of a single hot-spot also turned out to be irrelevant both for the initial and for long time epidemic behaviors.

In case of the 2d lattice model the diffusion had the effect of increasing the maximum of I(t) by ≃ 20% and doubling the duration with respect to the pure system at λ = 0.5. Further more detailed studies would be needed to clarify the hot-spot effects, involving a comparison with random sequential update model simulations.

It would be also be an interesting extension of research to clarify the effects of the intrinsic quenched disorder on the SIR dynamics as Harris criterion predicts irrelevance. Preliminary runs in d = 2 Euclidean lattices with 10% of randomly added quenched hot-spots do not modify the critical dynamics (see Appendix). However, the combination of hot-spots with mobility may affect the initial growth scaling more profoundly.

I thank Róbert Juhász for the discussions and Silvio Ferreira for the useful comments. Support from the Hungarian National Research, Development and Innovation Office NKFIH (K128989) is acknowledged. I thank access to the Hungarian National Supercomputer Network. pure case, with τ = 1.05(1) as can we see on Fig. 13 . Therefore, at least "weak" quenched interaction disorder seems to be irrelevant for the SIR critical behavior, not breaking the Harris criterion prediction. 

Nonequilibrium Phase Transitions in Lattice Models, Aléa-Saclay

Nonequilibrium phase transition: Absorbing Phase Transitions

Universality in nonequilibrium lattice systems: Theoretical foundations

Power-law decay in epidemics is likely due to interactions with the time-variant contact grap

Fractals and Disordered Systems

Chaos: An Interdisciplinary

In this Appendix I show some preliminary results for the study of interaction disorder in case of homogeneous, d = 2 Euclidean lattices of size L = 1000. The simulation methods are the same as discussed in Section III. I investigated heterogeneity in the form of quenched, intrinsic disorder. This was done by initializing the system with heterogeneous λ i -s. I used bi-modal disorder distribution