Diffusion-Localization Transition Point of Gravity Type Transport Model on Regular Ring Lattices and Bethe Lattices

Focusing on the diffusion-localization transition, we theoretically analyzed a nonlinear gravity-type transport model defined on networks called regular ring lattices, which have an intermediate structure between the complete graph and the simple ring. Exact eigenvalues were derived around steady states, and the values of the transition points were evaluated for the control parameter characterizing the nonlinearity. We also analyzed the case of the Bethe lattice (or Cayley tree) and found that the transition point is 1/2, which is the lowest value ever reported.

List of diffusion-localization transition points γ c for known synthetic graphs

Ring 1/2
Complete graph 1 Bidirectional star graph 1 Outward star graph ∞ Among these example graphs, the ring graph has the lowest γ c = 1/2 transport model is that the transport quantity depends on the power of the receiver's quantity, where the power exponent may be nonzero. Typical cases, including those related to thermal and chemical transport, do not exhibit this dependency, whereas nonzero power exponents are observed in some cases in the fields of social science and economics [10,14,15]. In this study, we focus on the effect of the power exponent on the transport properties. We define a parameter γ that controls transport nonlinearity. In a gravity-type transport system, larger quantities have a larger flow at positive γ values. Owing to this property, a point with a large quantity may monopolize a large portion of the flows for some γ . We call this phenomenon the diffusion-localization transition and focus on estimating the critical value γ c at which the stable steady-state solution bifurcates. Although some graphs with no feedback loops have a large γ c that exceeds 1, some graphs have a transition point γ c less than or equal to 1. A previous study [16] reported that a complete graph has γ c = 1, whereas a ring with infinitely many nodes has γ c = 1/2. Table 1 summarizes the γ c of several synthetic graphs. This paper introduces regular ring lattices as a solvable network model and derives an exact formula for γ c . We report the characterization of a structure with a low γ c in the context of the interaction length. We also consider the Bethe lattice as a typical graph to study mathematical models [17][18][19]. We analyze the system on the Bethe lattice both theoretically and numerically. The remainder of this paper is organized as follows. Section 2 describes the gravity interaction model used in this study. We also provide fundamental tools as the basis for our analysis in Sect. 2. In Sect. 3, we present an analysis of the diffusion-localization transition point on regular ring lattices as a generalization of the results reviewed in Sect. 2. In Sect. 4, we present the second main result of the analysis of the transition point on the Bethe lattices and present the master equation as an analog of that in Sect. 2. Finally, we discuss the results and clarify the correspondence between our analysis and the dynamical phase transition of a biased random walk in Sect. 5.

Model and Review of Examples
In [13], a nonlinear relation for Japanese firm-to-firm transactions is found such that the annual trading volume f i j from a buyer firm i to a seller firm j is proportional to the powers of their annual sales S i and S j , namely, f i j ∝ S α i S β j for some power exponents α > 0 and β ≥ 0; in addition, A = (A i j ) is the adjacency matrix of the inter-firm network such that A i j = 1 if there exists an edge from node i to j, and A i j = 0 otherwise. In the examples of world trade [10] and human flow between cities [11], gravity model with distance dependency has been empirically confirmed. However, unlike in the cases of world trade and human flow, in the case of firm-to-firm transactions [13], distance dependency is very weak, i.e., with δ nearly equal to 0, where d i j is the physical distance between i and j. From this relation, the following mathematical model on a directed network of N nodes is named as a generalized gravity interaction model [16]: where α and β are nonlinear power exponents, A = (A i j ) is the adjacency matrix of the given network, ν i S α i is the dissipation term that leaks out from the network, and F i is injected from outside the network. Here, t does not need to correspond to real time because we only need the stationary solution of Eq. (1), which is the realized S over the network. We assume that the network structure is fully described by its adjacency matrix A i j , whose value is 1 if there is a directed link from i to j and 0 otherwise. β = 0 corresponds to the equipartition where flow f i j does not depend on the sales of customer S j , also known as PageRank, and β = ∞ corresponds to a monopoly, where f i j is equal to S α i for j of the largest sales S j and 0 for other neighborhood nodes. In this paper, we only consider the case where F i is independent of both i and t, namely F i = F for some constant F. For the stationary solution S i of Eq.
(1) satisfies the following conservation law: We consider the following normalized quantity x i such that x i = S α i /F for the remainder of this paper. The equation for x i is where γ = β/α for α and β in Eq. (1). Hereafter, we set ν i = ν for all nodes i for simplicity. Thus, it is easy to confirm the following relation for the sum of the steady-state solution x from Eq. (2): In this model, the transport property is governed by the network structure A = (A i j ) and the nonlinear parameter γ . The formal statement of the diffusion-localization transition of this model is formulated as follows. There exists some γ c depending on the underlying network structure A = (A i j ) such that the steady-state solution of the model for γ < γ c bifurcates at γ = γ c . In other words, γ c is the value of γ at which the stable steady-state solution for γ < γ c becomes unstable. In this study, we used both theoretical tools and numerical simulations to study the transition point γ c . We start the numerical calculation from the uniform state 1/ν perturbed on one point to obtain a stationary solution when we perform numerical calculations. The perturbation will decay and diffuse over the system in the diffusion state, whereas the perturbation will grow in the localization state. For the diffusion state, we have the same stationary solution for different initial states, as far as we have performed the numerical calculation. For the numerical calculation, we determine the convergence by taking the L 2 -norm of the relative increment as less than or equal to 10 −6 . We have confirmed that the absolute error of conservation law (4) is less than or equal to 10 −2 under this convergence criterion.
Here, we review the previous results of the diffusion-localization transition on two illustrative examples of the ring and complete graphs. These two graphs have different transition points, 1 2 and 1. The diffusion-localization transition point γ c on the two graphs is derived by considering the master equation of the linearized equation.
We define the following N × N weight matrix of transport for k, i = 1, 2, · · · , N . The denominator is the sum over the out-neighborhood j of node k(see Fig. 1). We calculate the perturbation response of B based on the perturbation x from the stationary solution x. We consider cases in which a uniform steady-state solution is realized. When all node degrees are uniform, d, a uniform steady-state solution is realized by symmetry. In such cases, B i j = 1 d for adjacent i and j. For the uniform steady-state solution x i = x∀i over regular graphs, for a sufficiently small perturbation δx, the linear order of the perturbation response δ B ki of the matrix element B ki in Eq. (5) is for k, i = 1, 2, · · · , N . The last term represents the second-or higher-order terms of δx.
j represents the sum of the d neighborhood nodes of node k. Here, d is the coordination number; for example, d is N − 1 for the case of a complete graph, and 2 for the case of a ring.
Thus, from the model equation on the complete graph, the equation of perturbation is given by Thus, we can conclude that the transition point on the complete graph is 1 at N → ∞ for the following reasons: in the region of γ < 1, if δx i > δx j ∀ j = i, then dδx i dt < 0, which means that δx i decreases. On the other hand, in the region of γ > 1, if δx i > δx j ∀ j = i, then dδx i dt > 0, which means that δx i increases. This picture can be confirmed in Fig. 1b by observing the variance of x over the average of x scaled by size N . In the case of the ring, we calculate the equation of perturbation δx from a uniform steady-state solution as follows: We set r as the Euclidean coordinate of nodes i and r ± r as the coordinate of the neighborhood of node i. Then, we have the following set of equations as a result of the continuum limit: Here, we take the continuum limit of the equation to illustrate its relationship with the diffusion equation. From Eq. (9b), we find that the diffusion coefficient is positive for γ < 1/2 and negative for γ > 1/2. The first term, r 2 2 t , in Eq. (9b) is the usual diffusion coefficient in onedimensional space, and the second term, −γ r 2 t , is due to the nonlinearity of the interaction. This indicates that the diffusion-localization transition is caused by the two effects: diffusion in the usual sense determined by the network topology and the nonlinear effect of preferential flow towards nodes with larger x. Therefore, a phase transition between a normal diffusion phase and a localization phase occurs at γ c = 1/2. This diffusion equation establishes linearized gravity-type interactions as diffusion.
The diffusion-localization transition point γ c can also be calculated using a linear stability analysis. γ c in this sense is defined as γ , where the largest eigenvalue of the linearized matrix becomes positive. The general formula of the linearized matrix J = (J i j ) in Eq. (3) around x is given as follows: We omit the argument x of B to simplify the notation. The term governs the effect of nonlinearity on the linear stability. Note that the sum of the row entries satisfies j J i j = −ν for any x. Thus, J always has an eigenvalue −ν with the corresponding uniform eigenvector (1, 1, · · · , 1) ∈ R N . In the case of a complete graph, we have B i j = 1 N −1 for all i = j. Therefore, Thus, we obtain the exact eigenvalues of the linearized matrix given by Hence, γ c for finite N and ν can be calculated as follows: Thus, at a large size limit N → ∞ and zero dissipation limit ν → 0, γ c converges to 1.
On the other hand, in the case of the ring, we have B i j = 1 2 for j = i ± 1. Thus, the linearized matrix is given by Thus, the formula of the eigenvalues of circulant matrices [20] implies that for l = 0, 1, 2, · · · , N −1, depending on γ . γ c is the infimal value of γ such that max l λ l (γ ) > 0 by definition. Hence, γ c for the finite N and ν is As ν → 0, we can neglect the term ν/ sin 2 ( 2π j N ) in Eq. (14). This allows us to minimize j to 1.
Therefore, we have γ c → 1 2 at the large-size limit N → ∞. The Taylor expansion yields an asymptotic, which is described as where the symbol ∼ represents the major term in the limit of N → ∞. Eq. (16) was confirmed numerically, as shown in Fig. 1c. We confirmed that the diffused and localized solutions were stable. Figure 1d shows an example of diffused and localized solutions on the ring and their stability, starting from a uniform solution with perturbation on the origin. We obtained the uniform solution and localized solution on the ring for different γ = 0.000, 0.450, 0.518. In Fig. 1e, we show the decay of the relative change of the L 2 norm of the solution ||x(t + 1)||/||x(t)||. We further numerically confirm that the largest eigenvalue of the linearized matrix at γ > γ c is negative.

Analysis on Regular Ring Lattices
We consider regular ring lattices, as shown in Fig. 2a, as an example of the interpolation between the ring and the complete graph. The nodes are assumed to be arranged in a circle. Each node in this graph interacts with the 2k nearest neighbor. The model equation over a regular ring lattice is given by This equation includes the previous two examples, i.e., the ring for k = 1, and the complete graph for k = (N − 1)/2 as special cases. We have the exact eigenvalues of the linearized matrix on the uniform steady-state solution using the formula of the eigenvalues of the circulant matrices [20]. The linearized matrix of system (17) is derived as follows: The matrix element J i j depends only on the difference m = | j − i|. Therefore, we use both J m and J i j interchangeably.
Here, because m = | j − i| represents the hopping distance from i to j, 1 ≤ m ≤ k and k + 1 ≤ m ≤ 2k corresponding to distances from nodes i to j are 1 and 2, respectively. To derive Eq. (18), we count the 2-hop paths from nodes i to j to evaluate the third term of Eq.
depending on the value of γ . Note that l = 0 corresponds to λ 0 = −ν, which is the largest eigenvalue in the diffusion phase. From Eq. (19), γ c can be determined, as in the case of the ring. However, unlike in the case of the ring, we cannot theoretically determine l that maximizes λ l at γ = γ c . Thus, we employ a bisection method that seeks the infimal value of γ such that max l λ l (γ ) > 0, where the tolerance of the bisection method is set to 10 −6 . We obtained the transition point γ c , as shown in Fig. 2c. Figure 2c indicates that the transition point in the regular ring lattices continuously increases from ring (k = 1), where γ c = 1/2 to complete graph(k = (N − 1)/2), where γ c = 1 as k increases, and it shows an empirical relation where γ c is approximated by . This result reveals that the proportion of interacting partners within the total nodes determines the transition point. Eq. (19) includes the special case of Eq. (13) of k = 1.
The master equation is also calculated for a regular ring lattice. We use to calculate the following two equations: The master equation of perturbation (21) generalizes the two extreme cases of a complete graph and a ring. In the complete graph, the first and second terms of Eq. (21) can be combined such that Eq. (8) can then be recovered because 2k = N − 1 for a complete graph. The master equation of the ring can also be recovered from Eq. (21) by setting k = 1.

Analysis on Bethe Lattice
To determine the behavior of model (3) on the Bethe lattice, we assume "isotropic x" for all points in the same layer to ease the computational cost. This assumption enables us to map the Bethe lattice to the corresponding one-dimensional (1D) chain consisting of "radial layers" by polarization(see Fig. 3). Following the mapping, the r -th(r ≥ 1) layer of the graph consists of z(z − 1) r −1 equivalent nodes.
Using this mapping, we obtain the following set of equations from Eq. (3) for each point in the r -th layer of the mapped chain: When we perform a numerical simulation, we set the free boundary condition for the two outermost layers, that is, The steady states of Eqs. (22) and (23) satisfy the following conservation condition: As an approximation, we consider the following homogeneous equation without the dissipation and injection terms: Note that the uniform steady-state solution is one of the solutions to this equation. More importantly, however, the exponential solution x r ∝ (z − 1) −r is another solution to Eq. (25). Figure 4 shows an example of a steady-state solution for this equation. The steady-state solution of the model is almost uniform except in the neighborhood of the boundary layer for small γ , whereas an exponential-type steady-state solution is observed regardless of the system size for γ close to 1 2 . The uniform steady-state solution is approximately achieved for γ < γ c , except in the neighborhood of the boundary, while the exponential steady-state solution is approximately achieved near γ c . For the numerical calculation, we determine the convergence by taking the L 2 -norm of the relative increment as less than or equal to 10 −6 . We confirmed that the absolute error in Eq. (24) is less than or equal to 10 −2 under the convergence criterion.
We study the linearized equation of the model equation at a uniform steady-state solution with an infinite size limit. First, we sum up Eq. (22) for {x r } in the same layer r . Namely, for a set of the r -th layer points S r = {point in the r -th layer}, we define the total normalized quantity X r = i;i∈S r x i . Then, we obtain the following equation of transport between

Discussion
In this paper, we review an analysis of the diffusion-localization transition on the complete graph and the ring as prototypes and confirm that the transition point γ c is 1 and 0.5, respectively, at the large limit. Next, we generalize the results of the complete graph and the ring to the regular ring lattices by giving the general formula of the transition point γ c depending on the ratio of interacting pairs over the total number of nodes. We give a generalized eigenvalue formula of the linearized matrix, including the special cases of the complete graph and the ring, and show that regular ring lattices have transition points γ c between 1 2 and 1, corresponding to the ring and the complete graph, respectively. A challenge in generalizing the diffusion coefficient of the ring is that there is no available method for measuring the effect of the aggregation process in the neighborhood on the diffusiveness at a given site. This problem is essential because the structural dependency of the transition point γ c may be determined by such an effect. We further relate theoretical critical points 1 2 and 1 to the real-world example of firm-to-firm transactions reported in [16]. The transition point on firm-to-firm transactions is reported to be 0.9, between 1 2 and 1. This suggests that the link density of the transaction network is between the ring and the complete graph. We found the regular ring lattices as a series of graphs with a transition point between 0.5 and 1; in particular, we can construct a graph on which the transition point is 0.9. One further topic is clarifying the difference between this example and the real transaction network.
Finally, we show that the nonlinear gravity-type transport model on the Bethe lattice has γ c = 1 2 . By calculating the master equation and its drift coefficient of the linearized model, we argue that the diffusion-localization transition point γ c of the model at the large-size limit is 1 2 in networks other than Euclidean lattices. In our analysis, the key point is to map the system to a 1D chain by polarization. The diffusion-localization transition on the Bethe lattice is understood by a sign change in the drift coefficient of the master equation. This result verifies a microscopic picture of the diffusion-localization transition that particles at the origin that has gone outside once will return again. Unlike spin systems [18] in which the transition points of the 1D chain and Bethe lattice are different, their transition points of them are the same for the system we analyzed. Bethe lattice shows different transition points from the mean field results. The relation between the mean-field behavior and the Bethe lattice is a further discussion point to be studied. The methodology of our analysis using polarization is quite similar to that of the analysis of a biased random walk on the Bethe lattice [21], which assumes a constant uniform central field toward the origin. The model of biased random walk over the Bethe lattice with coordination number z is described as follows: the probability of hopping toward the origin is x z , where the bias parameter 0 ≤ x ≤ z and the probability of hopping to go further is z−x z . One difference between our model and the constant central field model is that our model does not directly assume a central bias toward the origin but assumes a localization parameter that conquers diffusiveness. The limitation of our approach is the assumption of polarization. The general formula of the transition point γ c on the regular ring lattices and Cayley trees based on linear stability analysis is still unknown. The behavior of the interaction length (see Fig. 2) and the asymptotic behavior of the transition point on the Bethe lattice (see Fig. 5) can be understood by finding such a formula to verify our numerical result in Fig. 2c of the regular ring lattices and the logarithmic convergence of the transition point of the Cayley trees to that of the Bethe lattice 1 2 . To study the diffusion-localization phenomena in real-world complex networks, we also need to consider the direction of the edges of the graphs. Because nonlinear gravity interactions have been observed in many fields, the mathematical foundation of this model, especially of the nonlinear effect of transport over a complex network topology, is expected to be studied.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.