The exact Laplacian spectrum for the Dyson hierarchical network

We consider the Dyson hierarchical graph , that is a weighted fully-connected graph, where the pattern of weights is ruled by the parameter σ ∈ (1/2, 1]. Exploiting the deterministic recursivity through which is built, we are able to derive explicitly the whole set of the eigenvalues and the eigenvectors for its Laplacian matrix. Given that the Laplacian operator is intrinsically implied in the analysis of dynamic processes (e.g., random walks) occurring on the graph, as well as in the investigation of the dynamical properties of connected structures themselves (e.g., vibrational structures and relaxation modes), this result allows addressing analytically a large class of problems. In particular, as examples of applications, we study the random walk and the continuous-time quantum walk embedded in , the relaxation times of a polymer whose structure is described by , and the community structure of in terms of modularity measures.

= − . L W J (1) Given that L is semi-definite positive, the Laplacian eigenvalues are all real and non-negative. Also, they are contained in the interval [0, min{N, 2w max }], where is called the Laplacian spectrum. Following the definition (1), the smallest eigenvalue ϕ 1 is always null and corresponds to the eigenvector μ 1 = e N , according to the Perron-Frobenius theorem. The second smallest eigenvalue is ϕ 2 ≥ 0, and it equals zero only if the graph is disconnected. Thus, the multiplicity of 0 as an eigenvalue of L corresponds to the number of components of .
The smallest non-zero Laplacian eigenvalue is often referred to as "spectral gap" (or also "algebraic connectivity") and it provides information on the effective bipartitioning of a graph. Basically, a graph with a "small" first non-trivial Laplacian eigenvalue has a relatively clean bisection (i.e., the smaller the spectral gap and the smaller the relative number of edges required to be cut away to generate a bipartition), conversely, a large spectral gap characterizes non-structured networks, with poor modular structure (see e.g., refs 30 and 31).
The spectral gap has also remarkable effects on dynamical processes. For instance, let us consider a system of N oscillators represented by the nodes of , which are interconnected pairwise by means of active links; the time evolution of the i-th oscillator is given by , where F and H are the evolution and the coupling functions, respectively, and β is a coupling constant. In this case a network exhibits good synchronizability if the eigenratio ϕ N /ϕ 2 is as small as possible and a small spectral gap is therefore likely to imply a poor synchronizability 32 .
However, probably the easiest dynamical process affected by the underlying topology is the random walk. In this context, the spectral gap is associated with spreading efficiency: random walks move around quickly and disseminate fluently on graphs with large spectral gap in the sense that they are very unlikely to stay long within a given subset of vertices unless its complementary subgraph is very small 33 .
Finally, we mention at the relation between the number of spanning trees in a graph  and its Laplacian spectrum. Let us denote with Ω( )  the number of spanning trees in : the Kirchhoff matrix-tree theorem states that (see e.g., refs 34 and 35) ϕ . This formula turns out to be extremely useful in order to get bounds for Ω( )  even without having the exact spectrum of the related graph, or in order to estimate the number of spanning trees of complex graphs which can be defined as combinations of simpler graphs for which the exact spectrum is known (see e.g., ref. 36).
For a weighted graph where  ∈ w ij is the number of edges joining i and j, the previous formula for Ω( )  still holds. More generally, when  ∈ w ij one can still look for the number of spanning trees on the bare topology (i.e., neglecting weights), or look for the minimum spanning tree(s), where the sum of the weights over all the edges making up the tree(s) is minimal (see e.g., ref. 37).
Scientific RepoRts | 7:39962 | DOI: 10.1038/srep39962 Beyond the purely mathematical point of view, deriving Ω( )  is a key problem in many areas of experimental design for the synthesis of reliable communication networks where the links of the network are subject to failure 38 . Also, the number of spanning trees is related to the configurational integral (or partition function) Z of a Gaussian macromolecule, whose architecture is described by ; in fact, one has ∝ Ω − Z [ ( )] 3/2  , where the proportionality constant does not depend on the topology but just on the number of monomers, on the temperature and on the spring constant between beads 39 .
Further applications for the Laplacian spectrum will be reviewed and deepened in the section "Examples of applications" focusing on the hierarchical Dyson network.

Description of the hierarchical network
In this work we focus on a deterministic, weighted, recursively grown graph, referred to as , originally introduced by Dyson to study the statistical-mechanics of spin systems beyond the mean-field scenario (corresponding to a fully-connected, unweighted embedding) 28 . The topological properties of this graph have been discussed in refs 14, 40 and 41, and here we briefly review them. The construction begins with 2 nodes, connected with a link carrying a weight J(1, 1, σ) = 4 −σ . We refer to this graph as  1 , and in the notation J(d, k, σ) we highlight the dependence on the graph iteration k and on the system parameter σ, also, d represents the iteration when the nodes considered first turn out to be connected. At the next step, one takes two replicas of  1 and connects the nodes pertaining to different replicas with links displaying a weight J(2, 2, σ) = 4 −2σ ; moreover, the weight on the existing links is updated as J(1, 1, σ) → J(1, 2, σ) = J(1, 1, σ) + J (2, 2, σ). This realizes the graph  2 , which counts overall 4 nodes. At the generic k-th iteration, one takes two replicas of − k 1  , insert 2 2k−1 new links, each carrying a weight J(k, k, σ) = 4 −kσ , among nodes pertaining to different replicas, and the weights on existing links are updated as If we stop the iterative procedure at the K-th iteration, the final graph K  counts N = 2 K nodes and the coupling between any pair of nodes can be expressed as Remarkably, this iterative procedure allows for a definition of metric: two nodes are said to be at distance d if they occur to be first connected at the d-th iteration [Note: One can check that this metric is intrinsically ultrametric since, beyond the standard conditions for a well defined metric (d ij ≥ 0; , the so-called ultrametric inequality (d ij ≤ max(d iz , d zj )) also holds.]. As a result, we can define a coupling matrix J associated to the network K  such that the entry J ij depends on the nodes (i, j) considered only through their distance d ij , namely Given the building procedure of  K , a generic node i has 2 d−1 nodes at distance d ∈ {1, … , K}. Moreover, the total weight of the links stemming from a single node i can be written as Due to the symmetry underlying the network, w i does not depend on the site i, so one can simply write w i = w. Beyond the coupling matrix J, one can introduce the Laplacian matrix L, which, as anticipated in Eq. 1, is defined as Before concluding this section it is worth stressing that the parameter σ is bounded as 1/2 < σ ≤ 1. In a statistical mechanics context this ensures that the Dyson model is thermodynamically well defined 13,14 ; in this context the lower bound σ > 1/2 ensures that the weight w remains finite in the limit N → ∞ , while the upper bound σ ≤ 1 ensures that the heterogeneity in the coupling pattern is not too strong (namely, that the relaxation time, given by the inverse of the spectral gap, does not grow faster than the system size).

Eigenvalues of the Dyson hierarchical graph
In order to get familiar with the structure of the matrix J, it is convenient to write it down explicitly for a small value of K. In particular, for K = 3 it reads as Scientific RepoRts | 7:39962 | DOI: 10.1038/srep39962 where we posed t = 4 −σ and we highlighted with the superscript (N) the size of the matrix; in general, N = 2 K and here N = 8. The block structure is evident and can be schematised as where L (N/2) and R (N/2) are square matrices of size × . For such matrices, the determinant can be computed as Using the ultrametric structure of J (N) , we can iterate this block decomposition. In fact, for example, after the first decomposition we obtain two matrices J N 1 /2 and J N 2 /2 such that ( /2) display the same block structure and for each an expression analogous to (7) holds in such a way that where, as stated previously, the superscript indicates the size of the matrices and, proceeding iteratively, we get a set of N/2 matrices of size 2 × 2 and referred to as More precisely, at every n-th iteration, we can write 2 n block matrices of size N/2 n , n = 1, … , K − 1, obtained as the sum or the difference of the 2 n−1 matrices of the previous iteration. With this argument, we can state that our determinant is a product of N/2 determinants of matrices 2 × 2 as Clearly, in order to compute the eigenvalues of J, we can again use this scheme, but taking as initial matrix J − λI. In this case, with some algebra, we obtain the following Scientific RepoRts | 7:39962 | DOI: 10.1038/srep39962 where p n is the determinant of a proper matrix of size N/2 n × N/2 n and it reads as By plugging the previous expression for p n into Eq. 11 and solving for the roots of det(J − λI) = 0, we get K + 1 distinct eigenvalues that can be written as K Notice that, as n increases, both λ n and its multiplicity decreases. Actually, we are mainly interested in the spectrum of the Laplacian matrix L = W − J, with W = wδ ij , where we recall that w is the weighted degree defined in Eq. 4. Now, denoting with ϕ n the n-th eigenvalue of L, and exploiting the fact that W is proportional to the identity matrix, one can write n n More precisely we have Of course, since L is semidefinite positive, ϕ n ≥ 0 ∀ ≥ n 0, and, in particular, omitting the trivial eigenvalue ϕ K = 0, we have From Eq. 16 we notice that the larger the value of ϕ n , and the larger its algebraic multiplicity (see also Fig. 1, panel a). Moreover, from Eq. 18 we see that ϕ 0 decreases with σ (see Fig. 1, panel b) and, as a consequence, the spectrum covers a narrower interval (we recall that ϕ K = 0 hence the spectrum width is given by ϕ 0 ). However, the eigenratio ϕ 0 /ϕ K−1 (which provides a more effective measure of the spectrum width) grows with σ implying poor synchronizability for large σ. This is consistent with the fact that large values of σ correspond to heterogeneous patterns of couplings 41 . Similarly, the spectral gap, corresponding to ϕ K−1 , decreases fast with σ (see Fig. 1, panel c) and, as a consequence, the spreading efficiency on K  is weakened as the patter of couplings gets more heterogeneous. Equations 16 and 17 provide an explicit expression for the spectrum of the Laplacian matrix associated to the hierarchical graph under study. The Laplacian spectral density ρ(ϕ), where ρ(ϕ)dϕ returns the fraction of eigenvalues laying in the interval (ϕ, ϕ + dϕ), can be analytically recovered as follows. First, we determine the cumulative distribution Cum(ϕ): fixing an arbitrary value ϕ ϕ ϕ ∈ [ , ] K 0 , the overall number of eigenvalues smaller than or equal to ϕ is given by where ϕ n( ) is obtained by inverting Eq. 16, that is, and where ⌊ ⌋ x denotes the largest integer not greater than x. We stress that ϕ ⌊ ⌋ n( ) provides the index n such that ϕ n is the eigenvalue closest to ϕ, with ϕ ϕ ≤ . Moreover, we remark that in Eq. 20, the index of the sum varies between ϕ ⌊ ⌋ n( ) and K − 1: this is because, according to the notation introduced in Eq. 16, the eigenvalues decrease as n increases, so to find all the eigenvalues smaller than or equal to a fixed ϕ, it is necessary to consider the larger indices.
Finally, taking the derivative of Cum(ϕ), we get an estimate for ρ(ϕ): n and, using (20) and (21), we get . This result is checked numerically in Fig. 3 (left panel): as expected, the comparison between the theoretical and the numerical estimates is successful especially in the case σ = 1, while for lower values deviations with respect to the power-law behavior emerge. Interestingly, as long as this picture holds, the expression in Eq. 24 allows us to get an estimate for the spectral dimension d s of . In fact, recalling that in the small eigenvalue limit ϕ → 0 one has ρ ϕ ϕ In particular, d s (σ = 1) ≈ 2, and it monotonically grows as σ approaches its lower bound. In fact, as σ is reduced, the pattern of couplings gets more and more homogenous towards a mean-field scenario 13,41 . This result is further checked in Fig. 3 (right panel). We stress that the spectral dimension (when defined) generalises the Euclidean dimension to the case of non-translationally-invariant structures: in fact, this is a global property of the graph, related (beyond to the density of small eigenvalues of the Laplacian), for instance, to the infrared singularity of the Gaussian process, or to the (graph average) of the long-time tail of the random walk return probability; it also provides a consistent criterion on whether a continuous symmetry breaks down at low temperature (see e.g., ref. 42).
Finally, it is worth comparing the results found here for  with those pertaining to other examples of graphs. For instance, if we neglect the weights in  we recover the complete graph (cg) of size N, often referred to as K N , where the Laplacian spectrum is given by ϕ (cg) = 0 with multiplicity 1 and ϕ (cg) = N with multiplicity N − 1. Beyond this peculiar case, as anticipated in the Introduction, there are several other examples of deterministically grown networks, exhibiting a non-trivial topology, for which the exact spectrum is known. In particular, we mention the Vicsek fractals (vf), where the spectral dimension is given by , with f being a parameter which sets the coordination number of inner points 43 , and a class of small-world networks (sw) obtained by recursively joining together K d graphs, where the spectral dimension is given by , regardless of d 26 . For relatively sparse graphs, such as exactly decimable fractals (e.g., Sieprinski gasket, T-fractal), one can prove that d s < 2 42,44 . The eigenvalue distribution ρ(ϕ)dϕ for a system of size N = 2 K with K = 10 is shown versus ϕ and for several choices of σ (the legend is the same as in Fig. 2). In particular, in this figure we show the histogram (symbols) derived from the exact spectrum, the theoretical estimate (solid line) given by Eq. 24, and the related scaling ρ(ϕ)dϕ ~ ϕ 1/(2σ−1) (dashed line). The comparison is successful especially for large values of ϕ and for σ not too small. Right panel: The eigenvalues for a system of size N = 2 K with K = 10 are plotted in ascending order, for several choices of σ, as reported. The dashed lines highlight the expected behaviour according to the analytical estimate ~(K − i) 2σ−1 , following from Eq. 22. Again, one can see that the theoretical picture is quantitatively good just for large values of σ, while it turns out to be qualitatively correct in the whole range considered. In fact, a slower rate of growth for ϕ i versus K − i implies a faster rate of growth for ρ(ϕ) versus ϕ, which in turns implies a larger spectral dimension, in agreement with Eq. 25.

Eigenvectors
Having computed the spectrum of L, we can now proceed with the calculation of its eigenvectors. As shown in the previous section, we have K + 1 distinct eigenvalues, each with its own algebraic multiplicity. We are looking for a system of linear independent vectors Again, we can exploit the block structure of L to find these eigenvectors. In fact, for K = 3, and posing again t = 4 −σ , L has the form Our ansatz is that is, a vector that has all the entries equal to zero, but those at the position 2j − 1 an 2j, that are respectively equal to +1/ 2 and −1/ 2 , in order to obtain N/2 vectors with norm equal to 1. In this case, the left hand side of Eq. 26 becomes Exploiting the structure of L, for every row i = 1, … , N, we have Comparing this expression with (28) we get that the ansatz (29) is correct. In this way we find a set of N/2 independent eigenvectors associated to the eigenvalue ϕ 0 .
We are going to proceed analogously for the computation of the eigenvectors related to for n = 1, … , K − 1, each with multiplicity N/2 n+1 . This implies that each of them is associated to 2 K−n−1 linear independent eigenvectors. In particular, we claim that they have the form Scientific RepoRts | 7:39962 | DOI: 10.1038/srep39962 where e is a one vector of size 2 n such that ||e|| = 2 n/2 and its first entry is either at the position 1 + i2 n , or at the position 1 With some algebra (as done previously in order to get (30) from (29)), we obtain where n = 1, … , K − 1, and i = 0, … , K − n − 1. The previous equation can be further simplified if one recasts the right-hand term of (32) as The last expression can be compared to ϕ j (see Eq. 16) times μ j (see Eq. 31), hence proving that the vectors defined in Eq. 31 are, in fact, eigenvectors. Finally, according to the Frobenius-Perron theorem, for the zero eigenvalue with algebraic multiplicity equal to one, the corresponding eigenvector In this way we have obtained a complete basis of N linearly independent eigenvectors related to the N eigenvalues, namely they form an orthonormal basis, that is i j i In particular, retaining the simple case K = 3, the orthonormal basis of eigenvectors can be written as The generalization to larger sizes is straightforward. We conclude this section observing that

Examples of applications
In this section we exploit the results found in the Sections "Eigenvalues of the Dyson hierarchical graph" and "Eigenvectors" to derive information about several processes which can be defined on the graph . These are just a few examples for illustrative purposes since, as underlined in the section "The Laplacian spectrum and its applications", the exact knowledge of the whole Laplacian spectrum can be applied in a very wide range of fields and situations.

Random Walks
A simple random walk embedded on an arbitrary graph  is characterized by the transition matrix P = W −1 J, that is, the probability that, in a given time step, the walker jumps from i to j is P ij = J ij /w i . Now, from the complete knowledge of the Laplacian spectrum one can derive the dynamical properties of the random walker and, in particular, one can calculate the mean time taken by a random walker to first reach a given node. This quantity plays a role in real situations such as transport in disordered media, neuron firing, spread of diseases and target search processes (see e.g., refs 44-46). Without loss of generality, we can fix the target site on the node j and focus on the mean first passage time (MFPT) from node i to node j, denoted as T ij . According to the definition of MFPT for random walks, we have T jj = 0 and, for any ≠ i j, 1 where e is the (N − 1)-dimensional unit vector (1, 1, … , 1) T ; T is the subvector of T obtained by deleting the j-th entry and whose i-th entry is T ij [Note: More precisely, the i-th entry of T is T ij as long as i < j, while if i > j the entry corresponding to T ij is the (i − 1)-th one]; P, W, and J are, respectively, the submatrices of P, W, and J obtained by deleting the j-th row and j-th column. With some passages it is possible to express T ij in terms of the spectra of L as (see ref. 27 for a detailed derivation) 1 . In particular, for the graph K  under study, we have that s = Nw, due to the homogeneity among nodes, and this also implies that T j is actually independent of j. Moreover, exploiting Eq. 36, we can further simplify the general expression (41) as The sum appearing in Eq. 43 is not feasible for an explicit, general solution, yet, noticing that  we can provide T j with a lower bound, that is, The reliability of this bound is shown in Fig. 4: by varying K and σ we see that the bound is more accurate when the size is large (i.e.,  K 1) and when the pattern of weights is less homogeneous (i.e., σ far from 1/2). The leading term of T lb scales as which evidences that (at least for large size and for σ far from its boundary values) the MFPT averaged over all starting nodes is well approximated by a linear growth with the size N, as corroborated by Fig. 5 (panel a). Moreover, one can show that (see also Fig. 4) j j and, for σ = 1, we are able to compute explicitly the sum in Eq. 42. In fact, fixing σ = 1, it becomes  where in the last approximation we outlined the leading term being C = 1/6 (see also Fig. 5, panel b). Therefore, by comparing the approximate result found for σ ∈ (1/2, 1) (see Eq. 46) and the exact result found for σ = 1 (see Eq. 49), we expect that T j depends functionally on σ.
In the case the location of the target is unknown, one can still obtain a characteristic time by averaging T j over all possible target locations, hence getting the so-called global first passage time, referred to as τ and defined as     where the last passage holds for large N (see Fig. 5, panel c).
To summarize, expected time T j (σ) to reach a target placed in the arbitrary site j on  can be bounded as j j lb

 
In general, as expected, T j (σ) grows with σ since, when the pattern of weights is more heterogeneous (i.e., σ is large), the graph tends to get disconnected, the mixing time is large and reaching the farthest nodes takes more and more time.
The leading asymptotic dependence on the network size displayed by T j (σ) can be compared to the scalings found for other structures displaying extensive connectivity (i.e., the coordination number of a subset of nodes diverges with the network size). For instance, for the class of pseudofractal networks studied in ref. 47, the MFPT Figure 5. Panel a: Behavior of log 2 (T j ) with respect to K = log 2 (N). Different sets of data, corresponding to different choices of σ, are shown with different symbols, as explained by the legend. Data points are obtained by numerical evaluating Eq. 43, while solid lines represent the best linear fit. In fact, as evidenced in Fig. 4, the lower bound given by T lb in Eq. 45 provides (in the limit of large size and for σ far from its boundary values) a good approximation for T j and, in turn, in the limit of large size, T lb scales linearly with N as shown in Eq. 46. More precisely, in our fitting function y = p 1 K + p 2 we let the fit coefficients free and in the inset we compare the best-fit coefficient p 1 with the unitary value expected from Eq. 46. The comparison is, within the error, satisfactory and further highlights that this approximation works better in the central region of the range σ ∈ (1/2, 1). Panel b: When σ = 1 an exact expression for T j is achievable and here we compare the asymptotical analytical result, reported in Eq. 49, with the numerical evaluation of T j from Eq. 43, as the size N = 2 K varies, being K ∈ [8,20]; notice the semilogarithmic scale plot. One can see that the numerical evaluation of T j (circles) is well approximated by the function CN log 2 (N), where C = 1/6 (solid line). Panel c: The expected time T j to reach a target set in the site j is plotted versus σ and for different choices of the graph size, as indicated by the legend. Data points are evaluated via the pseudo-inverse Laplacian method 70 (in order to get a further check), while solid lines represent the numerical evaluation of Eq. 43. Notice that, in the limit of large size, T j corresponds to the global first passage time τ.
scales sublinearly with the size N as long as the target node is central, while it scales linearly when the target node is peripheral; analogous results hold for the deterministic scale-free networks studied in ref. 48. Interestingly, a scaling ~N log N was evidenced for the global MFPT in the small-world exponential treelike networks studied in ref. 49, where a modular topology with loosely connections between different modules is also responsible for preventing the walker from exploring fluently the underlying space.

Quantum Walks
As stressed above, random walks constitute a basic example of dynamical process affected by the underlying topology. Their quantum-mechanical version, i.e. the quantum walks, constitute an advanced tool for building quantum algorithms and for modeling the coherent transport of mass, charge or energy (see e.g., refs 50-52), and display as well a strong dependence on the topological properties of the embedding structure and, accordingly, much of their properties can be expressed in terms of the Laplacian spectrum. Before deepening this point it is worth reviewing briefly a few definitions (here we just focus on the continuous time version of quantum walks, i.e., continuous-time quantum walks), while we refer to ref. 51 for an extensive treatment.
In such a quantum-mechanical context, the transport has to be formulated in Hilbert space and one assumes that the states |j〉 representing the nodes span the whole accessible Hilbert space; also, it is assumed that the states are orthonormal and complete, i.e., The dynamics is then governed by a specific Hamiltonian H, such that the Schrödinger equation for the transition amplitude α k,j (t) from state |j〉 to state |k〉 , reads where ħ has been set equal to 1. The squared magnitude of the transition amplitude provides the quantum mechanical transition probability π . Now, the quantum-mechanical Hamiltonian H can be identified with the classical transfer matrix, i.e., H ≡ L 53 .
The performance of the continuous-time quantum walk (CTQW) can be estimated in terms of the return probabilities π j,j (t): a quick temporal decay of these probabilities implies a fast transport through the network. In order to make a global statement on the performance, one considers the average return probability π given by (see ref. 51). where in the last passage we exploited the Cauchy-Schwarz inequality to get a lower bound for π t ( ) which does not depend on the eigenvectors.
For classical diffusion, the long-time limit of the transition probabilities reaches the equipartition value 1/N. In contrast, due to the unitary time evolution, for CTQW neither π i,j (t) nor π t ( ) decay to a given value at long times, but rather oscillate around the corresponding long-time average χ which, for π t ( ), is given by (see ref. 51) , lb n m namely, χ lb is the sum of the squares of the normalized multiplicities. Therefore, the full knowledge of the Laplacian eigenvalue spectrum allows estimating (at least via bounds) whether the quantum transport on a given structure propagates relatively fast. Notice that, actually, χ lb just depends on the multiplicities of the distinct eigenvalues. Indeed, for quantum transport processes the degeneracy of the eigenvalues plays an important role, as the differences between eigenvalues determine the temporal behaviour, while for classical transport the long time behaviour is dominated by the smallest eigenvalue. Situations where only a few, highly-degenerate eigenvalues are present are related to slow CTQW dynamics, while when all eigenvalues are non-degenerate the transport turns out to be efficient (i.e., it spreads out fast) 51 .
However, in general, an excitation does not stay forever in the system where it was created; the excitation either decays (radiatively or by exciton recombination) or, e.g., in the case of biological light-harvesting systems, it gets absorbed at the reaction center, where it is transformed into chemical energy. In such cases, the total probability to find the excitation within the network is not conserved. Such loss processes can be modeled phenomenologically by modifying the Hamiltonian H appearing in Eq. 54. To fix the ideas, we consider networks where the excitation can only vanish at certain nodes, making up a set  with cardinality = M  . These nodes will be called trap-nodes or traps. The new Hamiltonian reads as Γ ′ = − i H H , where Γ can be written as a diagonal matrix with elements Γ = Γ jj as long as ∈ i , and zero otherwise; Γ represents the (tunable) absorbing rate. As long as the absorbing rate Γ is small, i.e., Γ  1, this problem can be treated perturbatively getting that the mean survival probability Π t ( ) M can be approximated by a sum of exponentially decaying terms: Let us now resume the hierarchical graph K  and exploit the results obtained in Sections "Eigenvalues of the Dyson hierarchical graph" and "Eigenvectors" to derive an estimate for χ and for Π t ( ) M . We start with χ and notice that the estimate in Eq. 57 only accounts for the degeneracy of the eigenvalues which, as can be seen in Eqs 16 and 17, is not affected by σ, but only depends on the system size. As a result, χ lb does not depend on σ. Moreover, exploiting the expression for ϕ k obtained in Eqs 16 and 17, we can get that χ lb for a generic network of size N = 2 K reads as We can therefore derive that, no matter the system size and the value of σ, the long time average χ is always finite and larger than 1/3. This means that the coherent transport on K  tends to get stuck nearby the starting node. Of course, this stems from the fact that couplings decay exponentially fast with the distance among nodes. Actually, localization phenomena also occurs for the classical diffusion, where, as the size N gets larger and larger the mixing time grows exponentially and ergodicity breaks down in the limit N → ∞ 14,40,41 . However, at finite size, for classical propagation the equipartition is eventually reached, while for the coherent propagation the localization effect is significant already at small size.
Let us now move to Π t ( ) M ; according to Eq. 58 it is possible to get a long time estimate in terms of the imaginary part γ l of the eigenvalues of the perturbed Hamiltonian H′ . This can be calculated as follows: where in the second passage we exploited the diagonal form of the matrix Γ and in the third passage we exploited the results found in Section "Eigenvectors". More precisely, g is a vector, whose entry g l corresponds to the number of null entries in the eigenvector μ l , namely = … . − N g (0, 0, 1, 1, 2, 2, 2, 2, 3, , 2), and δ x,y is the Kronecker delta in such a way that the summation just returns the number of non-null entries in the eigenvector ψ l ( ) that match a trap position. For instance, when M = 1, one finds Moreover, the smallest eigenvalue is null and the second smallest eigenvalue is M/N, whose multiplicity grows with M. Each null eigenvalue corresponds to a block of K  which is trap free and this yields to a contribute to Π t ( ) M which survives even at long times. Therefore, the trap arrangement which maximizes the survival probability is the one where traps are set as close as possible each other, while the trap arrangement which minimizes the survival probability is the one where traps are scattered broadly. This is rather intuitive as the quantum walk tends to remain localized around its starting point in such a way that if traps are assembled within a block, the walker can avoid them as there exist eigenstates living on different blocks, while if each block (even the smallest one, i.e. the dimer) is occupied by at least one trap the quantum walk can not avoid them 56 .

Dynamics of polymer networks under external forces
Another field where the Laplacian spectrum is extensively exploited is polymer physics and, in particular, as for the investigation of the relationship between the topology of a polymer and its dynamics. To this purpose, the so-called generalized Gaussian structures are conveniently exploited. These are a generalization of the Rouse model 7,57 , meant for linear polymer chains, to systems with arbitrary topology. The polymer is modeled as a structure consisting of N beads (each bead corresponding to a node) connected by harmonic springs (we refer to ref. 7 for a discussion of the underlying assumptions).
Understanding how the underlying geometries of polymeric materials affect their dynamic behavior is becoming of increased importance as new polymeric materials with more and more complex architectures are synthesized 7 .
Let us consider a set of N beads, all subject to the same friction constant ζ with respect to the sorrounding viscous medium, and connected pairwise by harmonic strings. The pattern of the related string constants is provided by the matrix L in such a way that for the couple (i, j) it reads CL ij , where C is a suitable constant. The equation describing the motion of the j-th bead is the following Langevin equation where R j (t) = (X j (t), Y j (t), Z j (t)) represents the position in the space of the j-th bead at time t, L is the Laplacian matrix associated to the structure of the polymer, f j (t) is the thermal Gaussian noise [Note: More precisely, it is assumed that 〈 f j (t)〉 = 0 and ζδ δ δ where k B is the Boltzmann constant, T is the temperature, ζ is the friction constant, α and β represent the x, y, and z directions.], and F j (t) is the sum of all external forces acting on the j-th bead.
One can show that, acting (from time t = 0 onwards) with a constant external force F k (t) = Fδ jk on a single bead j of the polymer in the, say, y direction, the mean displacement of a bead turns out to be (see ref When the applied force is harmonic, namely say along the x direction, the related response function is the so-called complex dynamic (shear) modulus G*(ω), whose real, i.e. G′ (ω), and imaginary, i.e., G″ (ω), components are also referred to as the storage and the loss moduli, respectively 7 . In the generalized Gaussian structure model these are given by where ν represents the beads per unit volume.
Let us now resume the graph  K : as discussed previously, this represents a set of interconnected nodes where the coupling pattern exhibits hierarchy and modularity (see Fig. 6 Notice that the last term in the right hand side is proportional to the radius of gyration of the polymer considered 58 . The results collected in Eqs 66 and 67 show that, in the limit of short and long times the mean displacement varies linearly with time, with a rate which is independent of σ. Therefore, the most interesting regime is the intermediate one, where the pattern of weights can possibly play a role. In fact, as shown in Fig. 6, this is the case and the resulting displacement is enhanced when the network is more inhomogeneous (i.e., σ is large).
As for the storage and loss moduli, exploiting the results of Section "Eigenvalues of the Dyson hierarchical graph", Eqs 63 and 64 can be recast as  For several values of σ ∈ (1/2, 1], the behavior of G′ and of G″ has been evaluated numerically; the results for the case σ = 1 are shown in Fig. 7. As expected, G′ behaves as ω 2 at low frequencies and as ω 0 at high frequencies. In the intermediate regime the growth of G′ is smoothened and the width of such an intermediate scaling increases with σ. A similar behavior was evidenced for other unweighted networks displaying small-world features 26,59 . As for G″ , it peaks at around ω ≈ 1, the exact value depending on the choice of σ. For relatively small frequencies, G″ is larger when the pattern of weights is more inhomogeneous (i.e., σ large), while for relatively large frequencies G″ is larger when the pattern of weights is more homogeneous (i.e., σ small). In any case, the expected scalings G″ (ω) ~ ω and G″ (ω) ~ ω −1 are recovered. This is again consistent with the results found for small-world structures in refs 26 and 59.
For the particular case σ = 1 we can estimate the asymptotic behaviour of G′ by upper-bounding the expression in Eq. 68, that is The previous equation shows that G′ (ω)| σ=1 scales at most linearly with N. Moreover, in the limit of low frequencies (ω γ  N / ), from Eq. 70 we can get G′ (ω)| σ=1 ~ ω 2 N; this scaling is corroborated in Fig. 7 (left panel), where one can also see that this regime in the ω space shrinks with the size. In the opposite limit of high frequencies ω γ  ( 2 ), from Eq. 70 we get G′ (ω)| σ=1 ~ ω 0 N 0 ; again, the scaling is corroborated in Fig. 7 (left panel) and, this time, the width of this regime in the ω space does not depend on N.
Finally, as discussed in the section "Eigenvalues of the Dyson hierarchical graph", the case σ = 1 corresponds to d s = 2 which is a critical value over which the polymeric structures collapse (however, under external forces, such polymers can unfold 60 , which makes this analysis reasonable). For the case d s = 2 the observables studied in this section display characteristic asymptotic behaviors in the intermediate time/frequency domain (see ref. 26 and references therein), namely 〈 Y(t)〉 ~ log(t) and G′ (ω) ~ G″ (ω) ~ ω. These scalings have been successfully checked also for the graph  under study, as shown in Fig. 8.

Modularity
In graph theory the modularity Q is meant as a measure of the quality of a particular division of a network into a set of clusters (or groups, or communities) 29,61,62 . More precisely, for an unweighted graph, the modularity is the fraction of links falling within clusters minus the expected fraction in a network where links are placed at random (conditional on the given cluster memberships and the degrees of vertices). Thus, if the number of links within a group is no better than random, the modularity is zero. On the other hand, a modularity approaching one is indicative of a strong community structure, that is, a dense intra-group and a sparse inter-group connection pattern. Denoting with A the (unweighted) adjacency matrix of the graph, with ≡ ∑ z A i j ij the degree of a node i, and with i ij the number of edges in the graph, the modularity can be written as where  i is the cluster to which node i is assigned. The previous expression can be extended to assess community divisions in weighted networks as i ij is the overall weight spread over all the edges. Of course, Eq. 73 recovers Eq. 72 as long as weights collapse to the value 1.
Let us write Eq. 73 in a more convenient way 29 . First, we notice that the terms to be summed up in Eq. 73 can be considered as elements of the N × N matrix M given by  where η max is the largest eigenvalue of M. Now, the spectra of the modularity matrix is strongly related with the spectra of the weight matrix J (or, in the case of unweighted graphs, of the adjacency matrix A), see e.g., ref. 29. In particular, for regular graphs where = w w i , ∀i, the eigenvalues η { } i of the modularity matrix M are equal to the eigenvalues {λ i } of J (and, in turn, under a proper shift to eigenvalues {ϕ i } of L), but the largest eigenvalue λ K is replaced by the eigenvalue 0. This is proven rigorously for unweigthed graphs 29 and the proof can be straightforwardly extended to regular weighted graphs.
Specifically, for the graph  under study, recalling the results in Eqs 13 and 14, the spectrum of M reads as As examples of applications we studied the random walk moving isotropically in , the continuous-time quantum walk moving in a potential given by , and the relaxation times of a polymer whose structure is described by .
As for the simple random walk, we found that the expected time T j to first reach a given node j depends (super) linearly with the graph size N, the exact dependence being qualitatively affected by σ. In particular, when the pattern of weights is more inhomogeneous T j grows faster with N, as expected due to the unlikelihood to reach the farthest nodes.
As for the quantum walk, we proved that the long time average, measuring how spread the walker finally gets, is finite and lower bounded by 1/3 even in the limit of large size. This value has to be compared with the equipartition limit 1/N holding for the classical case, suggesting that the coherent transport is not well performing on . Also, in the presence of absorbing traps, according to their arrangement, surviving stationary states can be established.
Next, in the generalized Gaussian framework, we investigate the dynamics of polymers constituted by beads identified with the nodes of  and connected by elastic springs whose coupling is identified with the weights encoded by J. Using the Laplacian eigenvalues of  we estimated the mean displacement, the loss and the storage moduli under the stimulation of an external force. In particular, we found that in the intermediate time (or frequency) domain the response of the polymer is enhanced when the pattern of weights is more inhomogeneous. In the case σ = 1 we also recovered the expected scaling corresponding to a spectral dimension d s = 2.
Finally, we exploit our results on the Laplacian spectra in order to derive information about the modularity of the graph itself. In particular, we find that when the graph size is large the modularity is maximum (or close to its maximum value) for partitions of size scaling as the square root of the overall graph size.
The knowledge of the exact spectrum for  paves the way to further investigations in the phenomenologies related to systems embedded in hierarchical, modular structures. These analysis look particularly interesting given the peculiar behaviours evidenced, e.g., in cognitive and neuroscience [63][64][65] , and in bioinformatics [66][67][68] . For instance, the existence of a wide range of relaxation time scales may shed light on the existence of a number of metastable states for the ferromagnetic Dyson model 14 and on the feasibility of parallel retrieval in Dyson associative networks 69 .