Emergence of disconnected clusters in heterogeneous complex systems

Percolation theory dictates an intuitive picture depicting correlated regions in complex systems as densely connected clusters. While this picture might be adequate at small scales and apart from criticality, we show that highly correlated sites in complex systems can be inherently disconnected. This finding indicates a counter-intuitive organization of dynamical correlations, where functional similarity decouples from physical connectivity. We illustrate the phenomenon on the example of the disordered contact process (DCP) of infection spreading in heterogeneous systems. We apply numerical simulations and an asymptotically exact renormalization group technique (SDRG) in 1, 2 and 3 dimensional systems as well as in two-dimensional lattices with long-ranged interactions. We conclude that the critical dynamics is well captured by mostly one, highly correlated, but spatially disconnected cluster. Our findings indicate that at criticality the relevant, simultaneously infected sites typically do not directly interact with each other. Due to the similarity of the SDRG equations, our results hold also for the critical behavior of the disordered quantum Ising model, leading to quantum correlated, yet spatially disconnected, magnetic domains.

Scientific Reports | (2020) 10:21874 | https://doi.org/10.1038/s41598-020-78769-2 www.nature.com/scientificreports/ spreading model, sites with n i = 1 are referred as 'infected' , while sites with n i = 0 are 'healthy' or 'susceptible' . The contact process is a continuous-time Markov process on this state space, specified by the rates of possible (independent) transitions, which are the following (for an illustration see Fig. 1). First, infected sites become spontaneously healthy with a rate µ i , which may be site-dependent. Second, infected sites attempt to infect their adjacent sites with rates ij , and the trial is successful if the target site j was healthy. Again, the infection rates ij can vary from link to link. We will assume that the infection rate from site i to site j is the same as that from site j to site i, i.e. ij = ji . This variant of the model is also known as the susceptible-infected-susceptible (SIS) model 15 . This technical restriction, which is necessary for the applicability of the SDRG method, is irrelevant from the point of view of universal critical properties 16 .
The contact process exhibits a continuous, non-equilibrium phase transition between a non-fluctuating (absorbing) phase in which all sites are healthy and a fluctuating phase, where the density of infected sites is non-zero 17 . For regular lattices and uniform transition rates the transition falls into the robust universality class of directed percolation [18][19][20][21] .
If the transition rates are independent, random variables, the universality class is well characterized in one, two-and three-dimensional regular lattices due to large-scale Monte Carlo simulations [6][7][8][9][10] . The observed critical behavior is in line with the results of the strong-disorder renormalization group method [11][12][13][14] . According to this, the dynamical scaling of average quantities is ultra-slow, characterized by power laws of the logarithm of time rather than the time itself, whereas the static critical exponents are also altered compared to the directed percolation class 11,12,22 . Most interestingly, the critical exponents are independent of the form of the distribution of transition rates. Besides low dimensional regular lattices, a similar behavior has been found in spatially embedded networks with long-range connections [23][24][25] . These networks consist of a regular, d-dimensional, hypercubic lattice and a set of long links, which exist with a probability p ij ∼ βd(i, j) −s , where d(i, j) is the Euclidean distance between i and j. Note, that even for uniform transition rates these networks contain a so-called topological disorder due to the random connectivity of sites by long links. In the case s = 2d , the topological dimension is finite and varies with β [26][27][28] . As it was demonstrated in Refs. [23][24][25] for d = 1 and s = 2 , the contact process shows an ultra-slow scaling and the critical exponents vary with β , while an additional disorder in the transition rates is irrelevant. According to a general scaling argument presented in Ref. 29 , a qualitatively similar behavior is expected in higher dimensions for s = 2d.
Strong-disorder renormalization group. The strong-disorder renormalization group is the key method to understand the long-time behavior of the disordered contact process (DCP). For a general review and a detailed introduction we refer the reader to Ref. 13,14 . The first application of the method to the DCP was in Ref. 11,12 . Next, we recapitulate the essential features of the method.
The SDRG is a sequence of iterative steps operating on blocks of sites containing the largest rate (either an infection rate or a healing rate) of the model, see Fig. 2. If the largest rate is an infection rate, ij , the block of sites i and j is merged to a cluster, characterized by an effective healing rate: This simplification is a good approximation if µ i and µ j are much smaller than ij . If the largest parameter is a healing rate, µ i , site i is deleted, and effective infection rates are generated between all neighbors of site i: This approximation is justified under the condition ji , ik ≪ µ i . Apart from a one-dimensional lattice, it may happen that a newly formed cluster has, through its constituents, two infection rates to a third site. Similarly, if a new infection rate is generated, there may be a pre-existing infection rate between those two sites. These situations are usually treated by discarding the smaller infection rate, which is known as the maximum rule.
The above steps are applied sequentially, lowering thereby gradually the number of degrees of freedom, as well as the rate scale � = max{µ i , ij } . Regarding the set of clusters present in the system, the SDRG procedure can (1) lnμ ij = ln µ i + ln µ j − ln ij + ln 2.
(2) ln˜ jk = ln ji + ln ik − ln µ i . www.nature.com/scientificreports/ also be viewed as a special coagulation-annihilation process with extremal dynamics. The critical behavior of the DCP is governed by the so called infinite-randomness fixed-point (IRFP) of the SDRG transformation 11,12,30,31 , at least for sufficiently strong initial disorder 8,11,12,[32][33][34] . As the IRFP is approached along the critical renormalization trajectory, the distribution of logarithmic rates broadens without limits, and both types of reduction steps become asymptotically exact. Furthermore, this feature also justifies the applicability of the maximum rule, and even the neglection of the ln 2 term in Eq. (1). In the numerical SDRG calculations we therefore resorted to these approximations, which, as a further advantage, enable a computationally very efficient implementation of the SDRG procedure 35,36 .
Following the SDRG transformation of the DCP at a given realization of random rates down to a rate scale ≪ 1 provides an effective DCP with a smaller set of degrees of freedom, which approximates well the dynamics of the original system for times t ≫ −1 . Consider, for instance, the state of the system, evolved from a fully occupied initial state, at some late time t. Applying the SDRG procedure down to rate scale = t −1 , provides the set of sites with an O(1) occupation as the constituents of the clusters of the renormalized system, while the sites eliminated during the procedure will have an almost zero occupation probability.
In the case of a finite sample, the true steady state is always the absorbing, empty lattice state. Accordingly, all clusters are removed during the SDRG procedure, the last one at some rate 1 , the second last at � 2 > � 1 , in general, the nth cluster becomes irrelevant at � n > � n−1 . The inverse rates, � −1 1 , � −1 2 , . . . , are interpreted as the mean lifetimes of the corresponding clusters removed during the procedure. Clearly, the last one, −1 1 , gives the mean lifetime of the sample needed to be absorbed in the empty state. In finite but large systems in the vicinity of the critical point, the last few lifetimes such as −1 1 and −1 2 are typically well separated, meaning that they differ by orders of magnitude, and this time-scale separation becomes more and more pronounced with increasing system size. Consequently, in typical samples there is a considerable time span, −2 1 ≪ t ≪ −1 1 , within which practically no sites but the constituents of the lastly removed cluster are occupied.
The structure of the lastly removed cluster can also be used to determine a sample-dependent pseudo-critical point: making a double sized system by glueing together two copies of the original one, the onset of the active phase is indicated when the last cluster is different from that of the original system 22 .

Quasi-stationary simulation
The simulations were implemented for regular lattices in the following way. An occupied site (i) is randomly chosen and, with a probability µ/(µ + 1) it is made unoccupied, or, with a probability 1/(µ + 1) , one of its n neighbors (j) is picked with a uniform probability and infected with the probability ij , provided it was healthy. The time increment related to such a move is �t = 1/N s , where N s is the total number of infected sites.
For non-regular networks, where the sites may have different coordination numbers, a different implementation is needed in order to correctly simulate the process. Here, besides the list of infected sites, also a list of active links is stored, which contains all directed links with an infected source site. The number of elements of this list is denoted by N e . Then, with a probability µN s /(µN s + N e ) , a healing event occurs on one of the infected sites, chosen equiprobably, whereas, with the complementary probability, N e /(µN s + N e ) , an infection event is attempted. In this case, an active link (ij) is chosen with a uniform probability from the corresponding list, and the target site (j) is infected with a probability ij , provided it was healthy. The time increment of such an elementary move is �t = 1/(µN s + N e ).
We apply a quasi-stationary simulation with a reflecting boundary condition. This means that, at the point where there is only a single active site in the system, the healing event is rejected. This way we only need to validate that the system reached the quasi-stationary state, easily done by checking whether the mean density and its variance remain unchanged under increasing the relaxation time. The total simulation time was typically 10 28 , the first half of which is left for relaxation, while the measurements are performed in the second half.
First, in each sample, we determine an individual pseudo-critical point, which, for the system sizes we use, can significantly deviate from the ensemble average. Following the method proposed in Ref. 37 , we identify the pseudo-critical point with the maximum of susceptibility 38 defined as www.nature.com/scientificreports/ where N is the total number of sites, ρ denotes the global density of active sites, and �·� stands for the expected value in the quasi-stationary state.
Having estimated the pseudo-critical point, we performed here quasi-stationary simulations for a longer time, 10 30 , and measured the mean local densities (Figs. 3, 4, 5, 6). In addition, we started a number of independent simulations and averaged over them, in order to avoid the (rather improbable) possibility that the long-lasting activity is stuck at a cluster other than the one with the longest lifetime.
For a satisfying agreement with the SDRG method for moderate system sizes, the strength of the disorder needs to be sufficiently large, i.e. the distribution of the logarithmic infection rates needs to be sufficiently broad. In practice, the infection rates were chosen from a power law distribution, P < ( ) = 1/α , where the exponent α > 0 controls the strength of disorder.

Comparison of the cluster structure with the simulations
The average order parameter of the DCP is related to the fraction of original sites comprised by the clusters of the renormalized system, which decreases gradually as the SDRG proceeds. Its critical scaling properties can then be inferred from its dependence on the inverse time scale and inverse length scale (the number density Figure 3. Results in 1D systems. Asymptotic probability of being active in the quasi-stationary simulations as the function of the spatial coordinate for N = 1000 sites with periodic boundary conditions and strength of disorder α = 1 . As illustrated by the purple and red clusters, the SDRG is able to efficiently capture the activity profile even with a few clusters. The correlated clusters are disconnected fractals with a fractal dimension d f = 1+ √ 5 4 ≈ 0.819 30,31 , corresponding to a highly uneven activity profile. Gaps of strongly reduced density inside the clusters can be seen already at small scales, indicating asymptotically disconnected clusters of activity. Scientific Reports | (2020) 10:21874 | https://doi.org/10.1038/s41598-020-78769-2 www.nature.com/scientificreports/ of clusters) along the critical trajectory of the SDRG 13,14 . Instead of this, here we compare the spatial structure of SDRG clusters with the map of local occupancies obtained by simulations in individual samples. Ideally, the comparison of the two methods should be done under the same circumstances: (i) for the same random sample (set of random rates), (ii) at the same time (given by −1 in the SDRG method). The first requirement is, however, not the appropriate choice for a fair comparison of the methods. The reason is that the SDRG transformation, due to its approximative nature at early stages, does not preserve the position of the critical point. For instance, a truly critical initial system will depart from the critical trajectory and, in order to arrive at the critical IRFP, a compensatory initial shift in the control parameter is needed to be imposed. Therefore, rather than the initial point of the renormalization trajectory, its end point has to match the sample used in simulations. For this reason, the first requirement, i.e. the complete identity of rates must be relaxed, but a notion of "equivalence of the random environment" must still be preserved. Our approach to overcome this controversy was the following. We chose the infection rates randomly, while kept the healing rates constant. The latter can then serve as a control parameter and can be used to compensate the shift induced in the SDRG. The amount of the shift is quantitatively a priori unknown, therefore it is natural to choose the critical point as a "common point" to which the system can be tuned in both methods by using an indicator of criticality internal to that method. By this construction, the transition rates used in the two methods are not completely identical, but the difference is only in the uniform healing rate, while the random infection rates, which form the "random environment", and which govern the shaping of SDRG clusters, are identical.
Concerning the second requirement, we implemented the SDRG procedure up to the stage at which only one cluster remains in the system. Provided the lifetimes of clusters are well separated, this corresponds to the time span between the lifetime of the second last and the last cluster. The state of the system in this time span can be conveniently simulated, even without knowing the above lifetimes, by the above described quasi-stationary simulation, which prevents the system from getting absorbed in the empty state. www.nature.com/scientificreports/ In Fig. 3, we illustrate for three random samples that the main SDRG cluster (purple) not only contains all sites of nearly maximal probabilities, but also captures the majority of the probability weight. For such, moderately sized, systems smaller clusters still play a role, as shown by probability distribution captured by the secondary cluster. Interestingly, both the main and secondary clusters are fractured by deep gaps in the probability distribution, as signatures of asymptotically disconnected clusters. Overall, we see a good agreement already for moderate system sizes as illustrated in Figs. 3, 4, 5, 6 for d = 1, 2 and 3 dimensional lattices as well as for a two-dimensional lattice with long-range interactions (Fig. 6).

Discussion
In this paper, we studied the disordered contact process (DCP) in a random environment. By applying the combination of quasi-static simulations and an efficient renormalization group method, we showed that the critical behavior of the DCP is dominated by one strongly correlated cluster. Yet, as opposed to equilibrium systems, the governing clusters are predominantly disconnected objects, indicating that individuals who are infected at the same time typically do not know each other directly. According to the SDRG method, such disconnected clusters emerge from a strong, positive feedback mechanism, in which remote sites can effectively infect each other over and over again for a prolonged amount of time through indirect paths in the contact network 39 .
Besides infection spreading, variants of the DCP have recently gained increasing interest also in functional modeling of the brain 40 . Our results suggest that, as opposed to traditionally expected functional 'blobs' , strongly correlated brain regions are not necessarily directly connected. In other words, functional connectivity might qualitatively deviate from physical connectivity. This expectation exists in addition to the challenge that in the brain correlations do not necessarily decay with increasing physical distance 41 . We leave the extension of our results to dynamical models on experimentally obtained brain connectome datasets for future studies. Figure 5. Results in 3D systems. We again find a good agreement between the asymptotic density profile in the simulations and the clusters obtained by the SDRG method. In 3D, the fractal dimension is d f = 1.16(2) 35,36 , illustrated here for N = 50 × 50 × 50 sites with periodic boundary conditions and strength of disorder α = 4 . For better visibility, only sites with at least 10% of the maximum probability are indicated in the simulations, while we show only clusters of size above the square root of the size of the main cluster. Figure 6. Results in 2D with long-range interactions. Here we show our results for algebraically decaying interaction probability with an exponent s = 4 , for N = 100 × 100 sites with periodic boundary conditions. The number of short-cuts is N/16 and the strength of disorder is α = 4 . Red links indicate long-range connections within the main cluster. The critical clusters are even more strongly disconnected objects, with a formally zero fractal dimension 29 .