Dissecting localization phenomena of dynamical processes on networks

Localization phenomena permeate many branches of physics playing a fundamental role on dynamical processes evolving on heterogeneous networks. These localization analyses are frequently grounded, for example, on eigenvectors of adjacency or non-backtracking matrices which emerge in theories of dynamic processes near to an active to inactive transition. We advance in this problem gauging nodal activity to quantify the localization in dynamical processes on networks whether they are near to a transition or not. The method is generic and applicable to theory, stochastic simulations, and real data. We investigate spreading processes on a wide spectrum of networks, both analytically and numerically, showing that nodal activity can present complex patterns depending on the network structure. Using annealed networks we show that a localized state at the transition and an endemic phase just above it are not incompatible features of a spreading process. We also report that epidemic prevalence near to the transition is determined by the delocalized component of the network even when the analysis of the inverse participation ratio indicates a localized activity. Also, dynamical processes with distinct critical exponents can be described by the same localization pattern. Turning to quenched networks, a more complex picture, depending on the type of activation and on the range of degree exponent, is observed and discussed. Our work paves an important path for investigation of localized activity in spreading and other processes on networks.

Localization phenomena permeate many branches of physics playing a fundamental role on dynamical processes evolving on heterogeneous networks. These localization analyses are frequently grounded, for example, on eigenvectors of adjacency or non-backtracking matrices which emerge in theories of dynamic processes near to an active to inactive transition. We advance in this problem gauging nodal activity to quantify the localization in dynamical processes on networks whether they are near to a transition or not. The method is generic and applicable to theory, stochastic simulations, and real data. We investigate spreading processes on a wide spectrum of networks, both analytically and numerically, showing that nodal activity can present complex patterns depending on the network structure. Using annealed networks we show that a localized state at the transition and an endemic phase just above it are not incompatible features of a spreading process. We also report that epidemic prevalence near to the transition is determined by the delocalized component of the network even when the analysis of the inverse participation ratio indicates a localized activity. Also, dynamical processes with distinct critical exponents can be described by the same localization pattern. Turning to quenched networks, a more complex picture, depending on the type of activation and on the range of degree exponent, is observed and discussed. Our work paves an important path for investigation of localized activity in spreading and other processes on networks.

I. INTRODUCTION
Statistical mechanics of phase transitions was devised for infinitely large systems due to their singular behavior, possible only in the thermodynamic limit [1], in which subleading contributions from finite size regions to the order parameter become irrelevant. Dynamical processes that undergo nonequilibrium phase transitions share these same properties despite lacking a general background as in the equilibrium statistical mechanics [2]. However, the contribution of localized regions can change drastically the nature of the phase transition as, for example, its smearing which caused by extended disorder [3] in subextensive domains where the order parameter is locally nonzero but globally represents a negligible contribution. Another interesting phenomenon is the Griffiths phase, in which critical behavior in extended regions of the space parameter is produced by rare and locally active domains due to randomness in disordered systems [3].
Localization is central for activation transitions of dynamical processes on heterogeneous networks [4][5][6][7][8][9][10][11]. Let us consider dynamical systems where the agent states are represented by vertices of a network and their interactions by the edges connecting them. Examples include epidemic processes [12], opinion formation [13], and synchronization [14] models. One can define an order parameter that determines in which phase the system is. For example, a nonnull epidemic prevalence is the fraction of infected individuals that can indicate an endemic phase where the epidemic thrives in an extensive part of the system. Prevalence is an order parameter in the range [0, 1] that can be defined globally, locally, or even for a single vertex and represents an activity and can be generalized for other processes than epidemics. The local activity could be the probability that a vertex is infected in an epidemic or has the majority opinion in a voter-like model. Subgraphs such as hubs and their neighbors [15] or cliques of densely connected subgraphs [11] can hold high local activity for very long times (metastability) even if the overall order parameter is asymptotically null.
Dynamical processes on networks are frequently investigated within mean-field approaches [12,13,16], in which the transition between zero and nonzero order parameter is the main subject. When the transition is smooth, with the order parameter approaching zero continuously, these theories can be handled perturbatively at low prevalences using linear stability analyses [12,16] to investigate the behavior at the transition in terms of the spectral properties of matrices which are related to the network structure [4,[17][18][19][20][21]. A fundamental example is the susceptible-infected-susceptible (SIS) epidemic model, investigated thoroughly in the present work. In the SIS model, vertices can be infected or susceptible. Infected vertices become susceptible with rate µ and infect its susceptible nearest-neighbors with rate λ per contact. Let ρ i be the probability that vertex i = 1, . . . , N is infected and ρ = i ρ i /N the global prevalence. The evolution can be investigated using the quenched mean-field (QMF) theory encoded in the equation [4,12] in which A ij is the adjacency matrix defined as A ij = 1 if i and j are connected and A ij = 0 otherwise. The system undergoes a transition from an active to a disease-free state at the epidemic threshold where Λ 1 is the largest eigenvalue (LEV) corresponding to the principal eigenvector (PVE) v (1) = {v for λ − λ c λ c where β QMF = 1 and Since Λ 1 diverges as N → ∞ for random networks with power-law degree distribution of the form P (k) ∼ k −γ irrespective of the degree exponent γ > 2 [22], one has an asymptotically null QMF threshold for SIS on these networks. Based on Eqs. (3) and (4), Goltsev et al. [4] raised the question of whether the QMF prevalence near to the transition, when ρ i ∝ v (1) i , corresponds to an intensive quantity or not: An actual endemic phase demands that the PVE is delocalized and a 1 is nonzero whereas a localized active phase holds if a 1 → 0 as N → ∞. Note, however, that Eq. (3) is valid for λ − λ c λ c [4], representing a trick range since λ c → 0. In general, for finite networks the larger the spectral gap Λ 1 − Λ 2 , where Λ 2 is the second largest eigenvalue of A ij , the more accurate Eq. (3) is [17].
Goltsev et al. [4] used the inverse participation ratio (IPR) of the normalized PVE, defined as [4] to quantify localization in the QMF theory. The IPR assumes a finite value as N → ∞ for nonextensive localization on finite set of vertices, vanishes as Y 4 ∼ N −ν with 0 < ν < 1 for subextensive localization, and as Y 4 ∼ N −1 if the PVE is delocalized (extensive). According to this classification, the PVE for random networks is localized in a nonextensive region for degree exponent γ > 5/2 and in a subextensive region for γ < 5/2 [23], and in both cases corresponds to a localized phase. Moreover, the exponent β QMF = 1 in Eq. (3) is in odds with the exact result of Mountford et al. [24] where β > 1 for any γ > 2. However, recent numerical integration of the QMF equations and its extension to a pair approximation [25] show that these theories can fit simulations very accurately if not too close to the epidemic threshold, which is out of validity range of Eq. (3).
A message-passing theory for susceptible-infectedrecovered (SIR) model, for which an infected individual becomes recovered and cannot be reinfected (long-term immunization) [12], provides an epidemic threshold given by the inverse of LEV of the non-backtracking or Hashimoto matrix [26] constructed with the set of edges of the networks; see Refs. [20,26,27] for more details on Hashimoto matrices. In addition to the well-documented PVE localization of the adjacency matrix and its implication to the SIS dynamics, it was recently discussed localization in the PVE associated to the non-backtracking matrix on real networks and its implications for SIR and percolation processes [27].
The aforementioned results point out that localization analyses based on the spectral properties of particular matrices are limited and can lead to misleading conclusions with respect to the localization properties of the dynamical processes and the mean-field theories associated to them. Aspects such as dynamical correlations, stochastic fluctuations, and nonperturbative contributions cannot be fully reckoned in these spectral approaches. In the present work, we investigate localization of epidemic activity considering the actual dynamics on networks using a normalized activity vector (NAV) that gauges the contribution of each vertex of the network to the total activity. We analyze the SIS [12] and contact process (CP) [28] models due to their very distinct mechanisms of activation [29], but the methodology is not limited to them. Using both analytical and simulation tools, we show that the NAV components are directly associated to the vertex degree in annealed networks. Localization patterns for critical CP in annealed and quenched networks are essentially the same, corroborating that the critical behavior of this dynamical system is described by a heterogeneous mean-field (HMF) theory [28,30]. However, localization of critical SIS on quenched networks partially agrees with annealed networks for γ < 5/2 whereas highly structured and complex localization patterns emerge on quenched networks with γ > 5/2. We also show that CP and SIS on annealed networks present localization at the transition despite the delocalized, active phase that emerges just above it. Using fractional averages we show that the order parameter near to the transition is ruled by contributions of lowly localized vertices, which are by far the vast majority, while the IPR is determined by a vanishing fraction of the network. Finally, we analyzed the accuracy of the QMF theory in predicting the localization patterns observed in stochastic simulations of the SIS near to the epidemic threshold, and found that it is very accurate for γ < 5/2 but deviates quantitatively otherwise.
The remainder of the paper is organized as follows. Definition of the NAV, its determination, models, and other methods used in this work are presented in Sec. II. Effects of localization and their interpretation are exemplified applying the NAV tool to two simple networks models in Sec. III. The central core of the work, the analysis of localization on random networks with power-law degree distributions, is presented in Sec. IV. We finally draw our concluding remarks in Sec. V.

II. MODELS AND METHODS
Localization of dynamical processes on networks can be quantified using a local order parameter that measures the level of activity in a determined subset. We consider processes on the top of networks composed of i = 1, 2, . . . , N vertices or nodes and a set of edges or links that represent the interactions among them. The state of each vertex is represented by σ i and can be in one of two classes: active or inactive with probabilities ρ i and 1 − ρ i , respectively. We assume that the active vertices define the order parameter globally in the steady state (t → ∞) as which gives the stationary epidemic prevalence in both SIS and CP. The concept of active vertices is trivial in epidemic models with recurrent states such as the SIS. For the SIR model, the activity can be defined as the probability that a vertex will be infected in an outbreak started in a randomly selected vertex. For synchronization or opinion models the activity will depend explicitly on the neighborhood and should take into account higher order interactions such as pairwise. In this work, we deal only with the cases of activity defined at the nodal level since the generalization to higher order motifs can be done accordingly.
We define the activity vector ρ = (ρ 1 , . . . , ρ N ). We investigate processes with absorbing states on finite networks. We overcome the difficulties inherent to these systems by introducing a small self-activation [31] of vertices. It can be an exogenous infection in epidemics or, in general, spontaneous creation in reaction-diffusion processes. If the process has active steady states, self-activation can be introduced through a uniform and spontaneous rate f such that f → 0 as N → ∞ [32]. One advantage of this approach is that localization can be investigated also in the absorbing phase of the original dynamics. Obviously, other approaches, such as quasi-stationary simulations [32], can be used too.
In analogy to the PVE, we introduce the NAV φ = (φ i , . . . , φ N ) whose components are given by Two limit cases are the completely localized and delocalized states where ρ i = N δ ijρ and ρ i =ρ, respectively. The corresponding NAV components are φ i = δ ij and φ i = 1/ √ N , respectively. The IPR is defined as for which we derive the limit behaviors of Y 4 = N −1 for fully delocalized and Y 4 = 1 for fully localized activity. We analyzed localization in two recurrent activation models with binary states, namely, the SIS model [12] and the CP [2]. Contact process is a variation of the SIS model, in which the infection rate per contact is λ/k i where k i is the degree of the vertex which is transmitting the infection [28] while the healing process is exactly the same of the SIS. In both cases, states of a vertex can be represented by σ i = 0 when the vertex is inactive (susceptible) or σ i = 1 when it is active (infected). Despite similarities between the models, their activation processes are of very different nature when evolving on random power-law networks [29]. In the SIS, activation is triggered by hubs while in CP it happens collective as in an usual phase transition [1]. In the SIS on random networks with P (k) ∼ k −γ , hubs responsible for the activation can belong to a densely connected core for γ < 5/2 or be sparsely distributed for γ > 5/2 [15,23,29]. Thus, different patterns of localization are eligible for these models.
We investigated the models using HMF [33] and QMF [4] theories as well as stochastic simulations [34]. The stochastic simulations of the SIS with a source were performed using the following optimized Gillespie algorithm based on phantom processes [34]. At each time step, with probability one infected vertex is chosen at random and becomes susceptible. Here N inf is the total number of infected vertices and N SI is the total number of links emanating from them. With probability a vertex of the network is chosen with equal chance and becomes infected if it is susceptible. Finally, with probability 1 − Q − P an infected vertex i is selected with probability proportional to its degree k i and one of its neighbors j is chosen with equal chance. If j is susceptible it becomes infected. The time is incremented by where ξ is a random number uniformly distributed in the interval (0, 1). In CP simulations we replace N SI by N inf in Eqs. (9), (10), and (11) and chose with equal chance the infected vertex that transmits while all other steps are identical. See Ref. [34] for more details on phantom processes. We compute ρ in simulations as the fraction of time that each vertex i is active during an averaging time of t av after a relaxation time of t rlx . Mathematically it is expressed as Values t av = 10 5 µ −1 to 10 8 µ −1 and t rlx = 10 4 µ −1 to 10 6 µ −1 were considered, being the larger times used for lower activity regimes where fluctuations are larger. In numerical integrations, ρ is obtained after convergence to the steady-state. In all simulations we used f = µ/N meaning that at most one new infection is randomly introduced by one unit of time 1/µ.

III. LOCALIZATION PHENOMENA ON SIMPLE NETWORKS
In order to obtain insights about the localization in terms of the NAV, we start our analyses with simple networks A. Star graph A star graph consists of a central node i = 0 connected to i = 1, . . . , K neighbors of degree k = 1 (leaves). The adjacency matrix is A 0i = A i0 = 1 for i > 0 and A ij = 0 otherwise. Let us start with the QMF theory for SIS on the star graph defining ρ 0 as the probability that the center is infected and ρ 1 = ρ 2 = . . . = ρ N as the probability that a leaf is infected. The QMF equations become Without self-activation (f = 0) the epidemic threshold is λ c = µ/ √ K [35]. The steady-state solution for f µ and λ µ is given by and At λ = λ c = µ/ √ K we have ρ 0 = f /λ and ρ 1 = √ λf /µ leading to an IPR given by Y 4 = φ 4 0 + Kφ 4 1 , which asymptotically assumes Y 4 = 1/4 confirming localization. For supercritical region λ 2 K µ we have ρ 0 ≈ 1 and ρ 1 = λ/µ implying Y 4 = 1/K in a fully delocalized state. Finally, for subcritical region with λ 2 K µ 2 but still λK µ 1 in a subextensive localization since it goes to zero slower than K −1 . Figure 1(a) compares the IPR as a function of λ √ K/µ for the QMF theory of stochastic simulations of the SIS model on star graphs of different sizes K. The maximum localization in the QMF theory occurs at λ √ K/µ = 1 and the convergence to the asymptotic value Y 4 = 1/4 is verified. For simulations, the maximum value is Y 4 ≈ 0.08 occurring at λ √ K/µ ≈ 1.1 below the SIS threshold λ c √ K/µ = √ 2 for a pair QMF theory without self-activation [19]. In both cases, the localization drops very quickly after the transition since the range of λ where λ √ K/µ is finite corresponds to λ 1. The CP dynamics on star graphs has a finite lifespan for K → ∞ for any finite value of λµ [29] implying that the dynamics without a self-activation is always in the absorbing phase at long times. Developing a QMF theory for CP with self-activation, we have that ρ 0 = λC(µ, λ, f N ) and (17) is a finite positive constant when the arguments (µ, λ, ζ) are finite. Therefore, the IPR is asymptotically Y 4 = 1, as confirmed in Fig. 1(b) for both theory and simulations, showing that the CP dynamics with self-activation on a star is fully localized in the center for any finite value of (µ, λ).

B. Random regular networks with one outlier
Random regular networks (RRN) are simple graphs where all vertices have the same degree k i = m and connections are performed at random [35]. They are homogeneous networks where disorder is negligible. In this network, SIS and CP become equivalent by scaling the infection rate as λ SIS = mλ CP and we consider only the former model. The transition for f = 0 happens at a finite threshold, slightly above the pair approximation given by λ SIS c = 1/(m − 1) [19], and has a homogeneous mean-field like transition [2,36]. The SIS dynamics on a RRN is fully delocalized.
Localization can be introduced by the inclusion of a single vertex of degree K m to form a star subgraph immersed in the RRN [37]. The system can, in principle, present two activations being the first one of the star subgraph and other one of the remainder of the network at λ/µ ≈ 1/(m − 1). Figure 2 shows the epidemic prevalence and IPR as functions of the infection rate for a RRN with degree m = 6 plus a hub of degree K = 10 3 . We see two steep shoulders indicating the activations of the star subgraph and the RRN component. Remember that the transition is smoothed by the self-activation rate and would be singular only in the thermodynamic limit. Also observe the lower bound ρ = 1/N before the star activation. The IPR increases rapidly, reaching a maximum slightly before the first shoulder, and then decays slowly up to λ = λ RRN c , after which drops suddenly to a value Y 4 ∼ 1/N . Therefore, one can identify an extended localized phase in the range λ c < λ < λ RRN c , ending with a sudden delocalization. A very similar system, a Bethe lattice with a hub, was investigated in the framework of QMF theory [4], in which it was shown that the PEV is localized while the eigenvector of the second LEV is delocalized.
While insightful and providing the localization intensity, the IPR does not indicate where and how localization takes place. We can extract more information considering the NAV component probability distribution (NAVCPD) defined as the fraction of vertices Π(φ) whose components are in the interval φ i ∈ [φ, φ + δφ]. The NAVCPD for SIS dynamics near to the activation of the star subgraph is shown in Fig. 3(a). The distribution shows high localization in the hub (φ ≈ 0.35) and leaves connected to it (φ ≈ 0.029). However, the distribution also presents other spaced peaks representing, in principle, intermediary localization in successive shells of neighbors centered on the hub. We tested this hypothesis calculating the complementary cumulative probability (tail distribution) Θ(φ) defined as the probability that a vertex selected at random has NAV component larger or equal than φ. Figure 3(b) shows Θ(φ) corresponding to the bottom distribution of Fig. 3(a). We can clearly see the plateaus corresponding to the probability that a randomly selected vertex is a leaf, P leaf = K/N , their nearest, P nn ≈ (m − 1)K/N , and next-nearest, P nnn ≈ (m − 1) 2 K/N , neighbors. Observe that each plateau in Θ(φ) corresponds to a peak in Π(φ). Interestingly, other localizations are also detected in the NAVCPD, as indicated by other peaks and plateaus not associated with the aforementioned subgraphs. Finally, near to the activation of the RRN subgraph, top curve of Fig. 3(a), the NAV is concentrated around the delocal- ized value φ = 1/ √ N , but the distribution still presents localization patterns especially in leaves and center.

IV. LOCALIZATION ON POWER-LAW NETWORKS
We now consider synthetic uncorrelated networks with power-law degree distributions in the form P (k) ∝ k −γ with k = k 0 , . . . , k c . We consider both annealed and quenched networks with fixed minimal degree k 0 = 3. In the former, connections are rewired at random with a rate much higher than the typical rates of the dynamical processes evolving on them such that dynamical correlations are completely suppressed whereas in the latter edges are frozen and dynamical correlations are relevant [12].

A. Annealed networks
The classical HMF theory [33] becomes an exact prescription of annealed networks if fluctuations are negligible. The HMF equation for the probability ρ k that a vertex of degree k is infected is given by where Ω k is the probability that a randomly chosen neighbor is infected. We have Ω k = k P (k |k)ρ k for SIS [33] and Ω k = k P (k |k)ρ k /k for CP [38], in which P (k |k) is the probability that a vertex of degree k is connected to a vertex of k . The steady-state solution is For uncorrelated networks we have P (k |k) = k P (k )/ k [39] and Ω k = Ω independent of the degree. The epidemic threshold for SIS with f = 0 is given by λ c /µ = k / k 2 and goes to zero as size diverges only for γ < 3 [12] while for CP we have λ c = µ [40]. At λ = λ c we have λΩ ≈ ρ/ k = k P (k)ρ k / k in both SIS and CP models. This result can be obtained handling Eq. (19) in terms of hypergeometric functions and can be found in Ref. [33] for SIS and Ref. [40] for CP, for example. The NAV in terms of ρ k becomes while the corresponding IPR is The localization for λ = λ c is correlated with the vertex degree in HMF theory according to Eqs. (19) and (20): The higher the degree the more localized its activity is. According to the denominator of Eq. (19), the asymptotic behavior depends on k max ρ in the regime f ρ which we are considering. If k max ρ 1 the NAV becomes which leads to an IPR Equations (22) and (23) correspond to the PVE of the adjacency matrix of an annealed network given by A ij = k i k j /( k N ) presented, for example, in Refs. [4,23]. For example, the critical prevalence of the CP in the quasi-stationary (QS) state, where the dynamics returns to a previously visited configuration when the absorbing state is reached [32], goes to zero as ρ qs ∼ (gN ) −1/2 [41] where g = k 2 / k 2 . We numerically checked that the same result holds for CP and SIS on annealed networks using self-activation f = µ/N at the effective transition point given by the maximal of the dynamical susceptibility χ = N ( ρ 2 − ρ 2 )/ ρ [35]. For powerlaw degree distributions, we have k n ∼ k n+1−γ max for γ < n + 1 and k n ∼ const., otherwise, implying that ρk max ∼ N (γ−1−ω)/2ω for γ < 3 and ρk max ∼ N (2−ω)/2ω if γ > 3, in which it was assumed a general scaling for the maximal degree k max ∼ N 1/ω with ω > 1. We analyze natural ω = γ − 1 and structural ω = max(2, γ − 1) cutoffs [39,42]. The first one emerges when k c = N as a result of finite-size realizations of a power-law distribution while the second allows to generate uncorrelated networks for any γ > 2 [43].
Assuming ρ ∼ (gN ) −1/2 for a structural cutoff, we obtain k max ρ 1 for γ > 2 leading to implying that IPR vanishes for all values of γ > 2. Strong corrections to the pure power-law scaling are present for γ near to 3 and 5. Strict delocalization, with Y 4 ∼ N −1 , is found only for γ > 5 whereas subextensive localization, with Y 4 ∼ N −ν , and ν < 1 [23], is found otherwise. Interestingly, SIS and CP models have the same localization structure near to the transition whereas their critical behaviors for infinite networks with N = ∞ and k c = ∞ present different exponents for scale-free networks [33,40]. For the natural cutoff ω = γ − 1 and γ < 3, assuming again ρ ∼ (gN ) −1/2 , we have ρk max ∼ O(1) and in which F (a, b, c; x) is the Gauss hypergeometric function and n = 2 or 4. Taking the asymptotic limit with ρk 0 1 and keeping only the leading contribution in the first term we obtain in which Γ(x) is the Gamma function. Since ρk max ∼ O(1), so does F (n, γ − 1, γ; − k ρkmax ), implying that both terms in Eq. (26) are of the same order ρ γ−1 due to ρ ∼ 1/k max and ρ n k c n ρ γ−1 ∼ 1/N, (27) in which c n is a constant that depends on both epidemic model and degree distribution. Again, using Eq. (21), we have that the IPR becomes finite indicating nonextensive localization in a finite number of vertices. Note that while a subextensive fraction of the network is active due to ρ ∼ N −1/(γ−1) and, consequently, a diverging number of vertices is simultaneously infected as N increases, the localization identified by the IPR is only in a finite number of nodes. This issue will be discussed further in the sequence of paper. We performed simulations of both SIS and CP with selfinfection f = µ/N on annealed networks with power-law degree distributions using natural and structural cutoffs. Figure 4 shows the IPR analyses for γ = 2.3. Other values of γ are shown in comparisons with quenched networks; see Sec. IV B. The dependence of the IPR on infection rate is presented in Fig. 4(a). The CP presents a faster drop to zero at λ/λ c 1 in comparison with SIS model 1 . Notice that the abscissa is divided by λ c that goes to 0 and µ for SIS and CP, respectively. At the transition, Fig. 4(b), the predictions of the mean-field theory are confirmed. In the case of a structural cutoff, both SIS and CP present the same IPR value approaching zero according to Eq. (23). For γ > 5/2, the scaling given by Eq. (24) is confirmed but the prefactor is larger than the prediction of Eq. (23) as shown in Fig. 7. For the natural cutoff, the IPR becomes finite confirming the localization predicted by the mean-field theory. Finally, the NAVCPD presents peaks at each vertex degree present in the network as will be shown in the comparison with the quenched case; see Sec. IV B.
An interesting aspect calls our attention in this analysis. Activity in both CP and SIS dynamics in HMF or annealed approaches are localized at the transition according to the IPR analysis but an extensive, real endemic phase emerges above λ c [33,40]. So, we present examples where localization in the critical point does not discredit a meanfield theory for describing an endemic phase transition. A natural, but still unanswered question, is whether such a reasoning applies to the QMF theory as well.
A central issue for dynamic process with localization is to understand which subset of the network rules the epidemic prevalence and it was subject of intense studies [4,7,[44][45][46][47][48]. A misleading interpretation is that the epidemic prevalence is given by the most localized region indicated by the IPR as we are going to disentangle from now on. In other words, we are going to show that a finite IPR, which is dominated by a few nodes, does not mean that the epidemic prevalence is also dominated by these few nodes. Let the fractional average of a quantity F (φ) be defined as which yields the weight due to φ < φ in the average. For example, if F = 1 we obtain the cumulative probability distribution. Since φ i ∝ ρ i by definition, we have ζ{ρ n } = ζ{φ n } and ζ{Y 4 } = ζ{φ 4 }. The fractional average increases monotonically from 0 to 1 as φ varies from its minimal to maximal values. As an illustrative example, a fully delocalized state, in which Π = δ(φ − 1/ √ N ), has fractional averages given by a Heaviside function ζ = Θ(φ − 1/ √ N ). For HMF theory with a structural cutoff, the fractional average can be analytically computed using Eq. (22) to determine Π(φ) in terms of P (k) resulting and where φ 0 and φ max are given by Eq. (22) with k i = k 0 and k i = k max , respectively. Figure 5 presents the fractional averages computed for HMF theory with different degree exponents and k c = 2 √ N using Eqs. (29) and (30). The fractional averages for the IPRs become appreciable for φ near to φ max when ζ{Y 4 } ≈ 1, while for the prevalence it becomes appreciable just slightly above φ 0 and is close to 1 when ζ{Y 4 } is still departing from zero. The gap between fractional averages of prevalence and IPR increases with size implying that the epidemic prevalence is given by the delocalized component of the network. To show this, consider the value φ * for which the fractional average corresponds to a fraction (1 − c) of total average, in which c is small but finite (typically c = 0.1), yielding an estimate for the range of the leading contributions to the averages. The corresponding ratio of φ * for ρ and Y 4 for γ < 5 is This result should be understood as follows. While the IPR is dominated by a small, subextensive or nonextensive part of the system, relevant contributions for the order parameter come from the rest of the network. Moreover, as the network size increases the most active set contributes less for the overall activity and more for the IPR since k max diverges with N . The interpretation of localization near to a transition is made precise: While a small portion of the network can have extremely high activity that rules the IPR analysis, its contribution to the epidemic prevalence is negligible; the latter is given by the rest of the network that has lower level of activity but is much more numerous.

B. Quenched networks
The simpler analytical toolbox based on HMF theory used for annealed networks is not available for the quenched case. Therefore, we use stochastic simulation [34] and numerical integration of the QMF equations [25]. For the latter, we include the self-activation term in Eq. (1) to obtain (32) which was solved numerically using a fourth order Runge-Kutta method. The QMF equation for CP can be obtained from Eq. (32) replacing A ij by A ij /k j . We used the uncorrelated configuration model (UCM) [43] with a structural cutoff k c = 2 √ N to generate quenched networks. Figure 6 shows the NAVCPD for simulations of the SIS and CP dynamics at their corresponding transition points running on UCM networks with γ = 2.3. Both models present very similar values of the IPR, consistent with the decays for annealed networks given by Eq. (23); see inset of Fig. 4(b). While the distribution for CP presents peaks corresponding to φ i ∝ k i and matching almost exactly the annealed counterpart, the distribution for SIS matches almost perfectly the QMF theory, agreeing with the annealed case only in the tails. The QMF theory for CP at the transition also presents the peaks φ i ∝ k i . So, the IPR alone fails in distinguishing differences in the localization of these dynamic processes. As discussed in Section II, SIS and CP on quenched networks are characterized by different activation mechanisms [29,49] and differences in the localization are expected.
The NAVCPD for critical SIS on quenched networks with γ = 2.3 presents two asymptotic regimes: It is approximately Π ∼ φ −1 for low localization while the tail scales as does the degree distribution Π ∼ φ −2.3 . We numerically determined that the crossover point between these regimes is very close to in which the equality can be obtained from Eq. (7). This crossover is related to the nature of the activation of the SIS process for this range of γ, which is triggered in a densely connected component of the network identified by the maximum index of a k-core decomposition [15]. The k-core decomposition is the removal of vertices and edges connected to them such that only vertices of degree k ≥ q remain in the network, where q is the k-core index [50]. The maximum k-core corresponds to the last step before all vertices are removed. Figure 6(a) shows the NAVCPD computed for SIS only on a subgraph containing the vertices belonging to the maximum k-core plus their nearest neighbors. The distribution in this subgraph matches very well the overall distribution for φ > φ * , dropping quickly for φ < φ * , which corroborates our assertion. Note that the maximal k-core itself corresponds only to the narrow end of the tail. This subset corresponds to a vanishing fraction of the network as N → ∞ [49] that dominates the contributions for IPR; see discussion on Fig. 8. Some important consequences come from Fig. 6 and analyses of CP for other values of the degree exponents. A first one is the corroboration that CP on quenched networks is a HMF-like process [28,30], an issue intensely debated in the 2010s [38,40,[51][52][53], due to the almost perfect match of NAVCPD for simulations on quenched and annealed networks. For SIS, a remarkable agreement between NAVCPD for QMF theory and simulations occurs showing that the QMF is able to capture almost perfectly the localization pattern observed in simulations for the investigated range size (up to N = 10 7 ), even though the critical exponents of the prevalence are different with β > β QMF = 1 [21,24]. Notice that a similar situation happens for CP and SIS on annealed networks as discussed in Sec. IV A. Obviously, one cannot discard the possibility of a different scenario for much larger, computationally inaccessible sizes. Finally, the crossover in the NAVCPD of SIS also indicates that the HMF behavior holds for an extended maximum k-core that encloses its nearestneighbors.
For γ > 5/2 striking divergences between SIS dynamics on annealed and quenched networks emerge [15,24,35], the more evident for the larger γ. Figure 7 shows the NAVCPD analyses of the critical SIS for γ = 3.5 on both quenched and annealed networks with a same degree distribution. For annealed networks, the behavior is similar to the case γ = 2.3 with peaks proportionally related to the vertex degree. The corresponding IPR decays according to Eq. (23) as Y 4 ∼ N −0.4 , but the prefactor is not accurate as it was for γ = 2.3; see Fig. 6(b). Simulations on annealed networks are less localized than the HMF theory?s prediction due to the fluctuations absent in theory which contribute to spread the activity reducing, therefore, the localization. In the case of quenched UCM networks, the IPR seems to converge to a finite value but it has not achieved its asymptotic value in the investigated size range with Y 4 ≈ 0.01 for N = 10 7 . The NAVCPD for quenched networks differs substantially from their annealed counterparts. No correlation between NAV and degree distribution is evident. Also, it does not quantitatively match with QMF theory either, but it seems to capture the trend; see discussion on Fig. 9. The NAVCPD for the quenched network presents very heavy tails, decaying approximately as Π ∼ φ −η with η 2 (the actual η value is beyond our scope), being therefore ruled by outliers that determine the IPR value. The case γ = 2.8, which is not shown for brevity, is similar to γ = 3.5. The fractional averages of prevalence and IPR for SIS at λ = λ c are presented in Fig. 8 for single realizations of UCM networks with N = 10 7 vertices and three values of the degree exponent. As in the annealed case, we have that vertices of lowest φ values contribute significantly for prevalence while the IPR is determined by the most active ones; the latter contribute negligibly for the prevalence value. Curves for γ = 2.3 are smooth and qualitatively similar to the annealed case, with significant differences only for prevalence at very low φ. In turn, curves for larger values of γ contrast with the annealed case, being more evident for larger γ, in which the fractional IPRs vary abruptly for large φ. In Fig. 8, for example, the vertex of highest φ represents 22% and 66% of the total IPR for networks with γ = 2.8 and 3.5, respectively. Curve ζ{ρ} for γ = 3.5 presents several jumps and plateaus. The same happens for the corresponding IPR, which can be seen in a (not shown) double-logarithm scale. This is an explicit observation of localized activity that is known to drive the epidemic activation for SIS at UCM networks with γ > 3 [29,45,54].
Important insights and elucidations are extracted from Figs. 6, 7, and 8, helping to rationalize whether metastability of weakly interacting subgraphs leads or not to an actual endemic phase [4,11,[44][45][46][47][48]. The SIS dynamics on quenched networks at the epidemic threshold is indeed localized for all values of γ (subextensively for γ < 5/2 and on a finite set for γ > 5/2) as does the PVE of the adjacency matrix [23], in the sense that a vanishing minority of vertices presents extremely high epidemic activity, ruling the IPR analysis. However, the corresponding epidemic process can still be undergoing a phase transition, in which extensive components are participating in the processes. Figure 9(a) presents log-binned scatter plots for NAV components obtained in SIS simulations (φ) and integration of the QMF equations (φ QMF ) for an infection rate λ = λ c determined by the maximum of χ. In contrast with the case γ = 2.3, for which an almost perfect match is observed (see Fig. 6), the scatter plots for γ = 3.5 present a very disperse distribution of points. The integration captures the trend of the simulation presenting, however, most vertices with NAV components higher than in stochastic simulations, which is indicated by most points below the identity line. This is due to the QMF equation being actually supercritical since λ c > λ QMF c . In Fig. 9b), we compare the situation in which both QMF theory and stochastic simulations are in their corresponding epidemic thresholds (λ = λ c for simulations and λ = λ QMF c = 1/Λ 1 for QMF theory). Here, the converse of Fig. 9(a) is observed: Activity in QMF theory is slightly correlated with stochastic simulations only for the few and most active vertices whereas the vast majority of the network presents no correlation.
In order to illustrate the distinct levels of localization in these approaches, consider the single vertex contributions to the IPR in the network sample shown in Fig. 9. The vertex indicated by arrows is the one of higher activity in both simulations and QMF theory. However, the gap between the first, φ 1 , and second, φ 2 , most active vertices is much larger in critical QMF than in simulations, given by φ 1 /φ 2 ≈ 15 and 1.8, respectively. In both cases, the IPR is led by the most localized vertex but much more in the QMF theory. Indeed, this single vertex contributes with (φ QMF 1 ) 4 = 0.2380 to the IPR Y 4 (φ QMF ) = 0.2384 for critical QMF theory, and with φ 4 1 = 0.0074 of Y 4 (φ) = 0.0087 for stochastic simulations. For supercritical QMF analyzed in Fig. 9a), we have (φ QMF 1 ) 4 = 0.0007 of Y 4 (φ QMF ) = 0.0029. Note that the IPR of critical QMF is very close to the value Y 4 = 1/4 of a star graph; see Section III A.
From the perspective of theory accuracy, these results show that QMF theory should be used with extreme caution to determine the most active vertices in the network since it may deviate substantially from the actual picture as shown in Fig. 9a). An alternative is to use pair QMF theory [19,25] which reckons dynamical correlations considering pairwise equations. Indeed, a pair QMF theory was recently used to identify the most central spreaders that should be immunized to efficiently stop a SIS dynamics [55]. However, even pair QMF theory has limitations for SIS dynamics on networks with large γ [21] and should also be used with caution.

V. CONCLUSIONS
Effects of localization introduced by structural disorder play a major rule on equilibrium and nonequilibrium statistical mechanics [3] and have recently been applied to dynamical processes on networks [4-6, 9, 10], in which disorder is intrinsic to the heterogeneous connectivity structure observed in real and synthetic networks. The activity localization in dynamical processes has been widely, not exclusively, investigated in terms of spectral properties of Jacobian matrices that emerge in the linearization of mean-field equations near to the transition [4,17,21,23,56,57]. In the present work, we develop a localization analysis applicable to any type of dynamical process for which activity can be gauged in terms of a local order parameter which, in present work, was considered at a vertex level but lower resolution motifs are eligible as well. The method is generic and can be applied to both mean-field and stochastic simulations whether the system is near to a transition or not.
We applied the methodology to the SIS [12] and CP [38] models, two conceptually similar dynamical processes with very distinct activation mechanisms, to a wide spectrum of random networks which were handled both analytically (mean-field theories) or using simulation tools. We were able to disentangle some misleading interpretations of localization in epidemic models. One conclusion is that the IPR, a commonly used metric to quantify localization in an overall level, blurs the structured localization patterns present in complex and even relatively simple networks.
Handling the HMF theories for SIS and CP on networks with power-law degree distributions, we have highlighted several elucidative outcomes. We show that SIS and CP present exactly the same localization patterns at the epidemic threshold for any value of the degree exponent γ, despite these models being characterized by different critical exponents for the epidemic prevalence when γ < 4 [33,38].
Also, we show that both CP and SIS dynamics at their respective thresholds are localized in a subextensive subset (increasing sublinearly with network size) for γ < 5 and concomitantly give rise to an endemic (delocalized) phase immediately above the threshold. Finally, defining the concept of fractional averages we show that even in states which are localized according to the IPR analysis, the order parameter is determined by the vast majority of vertices in which the activity is smaller in contrast with the IPR which is ruled by a fewer (subextensive), small fraction of the network. Even seeming unsound at a first glance, the interpretation is quite simple. A few vertices are active most of the time while others are active only eventually due to stochastic interactions with the former. However, the latter being infinitely larger in number will lead any intensive quantity such as the order parameter. These finds are foundational for the correct interpretation of activation mechanisms of endemic phases.
We also analyzed localization of both dynamics on quenched networks, which are in general not accurately described by mean-field approaches and were handled using stochastic simulations. Firstly, we corroborate the HMF-like nature of the CP on quenched networks, an issue widely debated formerly [38,40,[51][52][53] in terms of meanfield exponents, showing that annealed and quenched versions of a same degree distribution have the same localization patterns depending only on the vertex degree. We stress again that the localization pattern alone cannot give the ultimate response and it complements other evidence reported elsewhere [28,30].
The results for SIS on quenched networks are richer and more complex. For small γ < 5/2, the localization patterns at the threshold match almost perfectly the QMF theory, agreeing with the annealed case only for a subset essentially composed by the vertices of the maximum kcore plus its nearest neighbors, being still a subextensive subgraph. Again, one knows that critical exponent of the epidemic prevalence of the QMF theory β QMF = 1 is not the correct one since β > 1 for any value of γ > 2 [24]. Despite this, our numerical results suggest that QMF is still able to capture accurately the localization structure. For γ > 5/2, the distribution of activity localization is more complex presenting outliers and heavy tailed distributions that lead to epidemic localization in a finite subset, as predicted by the QMF theory without however an accurate quantitative match with the simulations. Indeed, comparing critical behavior of both simulations and QMF theory, correlation is found only for a tiny part of the vertices in which activity is highly localized.
To sum up, our results show that the localization patterns are much more complex and revealing than simple metrics such as the IPR can tell. Using spectral properties of, for example, adjacency or non-backtracking matrices, derived within mean-field frameworks, can blur important features of the dynamical process and also lead to misleading or incomplete conclusions on the physical mechanisms involved in. The methodology studied here can be extended to any class of networks and dynamical processes. Therefore, we expect many forthcoming analyses of localization of other dynamical processes, such as opinion dynamics [13] and synchronization [14], as well as other network topologies such as multilayer [8] and temporal [58] networks.