Diffusion dynamics and synchronizability of hierarchical products of networks

TL;DR

This paper investigates how hierarchical network structures influence diffusion and synchronization, revealing distinct regimes and optimal coupling conditions that depend on the subnetworks' properties.

Contribution

It introduces a method to analyze eigenvalues of the Laplacian for hierarchical networks, elucidating how primary and secondary subnetworks affect dynamics.

Findings

Diffusion rate depends on secondary network at small coupling.

Diffusion saturates at large coupling, dominated by primary network.

Synchronization is optimized at an intermediate coupling value.

Abstract

The hierarchical product of networks represents a natural tool for building large networks out of two smaller subnetworks: a primary subnetwork and a secondary subnetwork. Here we study the dynamics of diffusion and synchronization processes on hierarchical products. We apply techniques previously used for approximating the eigenvalues of the adjacency matrix to the Laplacian matrix, allowing us to quantify the effects that the primary and secondary subnetworks have on diffusion and synchronization in terms of a coupling parameter that weighs the secondary subnetwork relative to the primary subnetwork. Diffusion processes are separated into two regimes: for small coupling the diffusion rate is determined by the structure of the secondary network, scaling with the coupling parameter, while for large coupling it is determined by the primary network and saturates. Synchronization, on the other hand, is separated into three regimes: for both small and large coupling hierarchical products have poorly synchronization properties, but is optimized at an intermediate value. Moreover, the critical coupling value that optimizes synchronization is shaped by the relative connectivities of the primary and secondary subnetworks, compensating for significant differences between the two subnetworks.

Equations

$$A_{\alpha }$$

$$L_{\alpha }$$

$$L_{\alpha }({\mu }_i)=L_1+\alpha {\mu }_i{D}_1,$$

$$L_{\alpha }({\mu }_1=0)=L_1.$$

$$\lambda (\alpha )={\nu }_j+\alpha {\mu }_i{\mathbf{v}}^{jT}{D}_1{\mathbf{v}}^j,$$

$$\lambda (\alpha )=\alpha {\mu }_i+{\nu }_j^{\cancel{0}},$$

$$\lambda (\alpha )={\nu }_j^0,$$

$${\dot{x}}_i=\sum _{j=1}^N{a}_{ij}(x_j-x_i),$$

$$\dot{\mathbf{x}}=-L\mathbf{x},$$

$$\displaystyle\|\bm{x}(t)-\bm{x}_{\infty}\|\mathrel{\vbox{ \offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern 2.0pt\cr\sim\cr\kern-2.0pt\cr}}}e^{-\lambda_{2}t},$$

$${\lambda }_2(\alpha )=\alpha {\mu }_2{\mathbf{v}}^{1T}{D}_1{\mathbf{v}}^1=\frac{\alpha {\mu }_{2}n}{N_1},$$

$${\lambda }_2(\alpha )=\text{min}({\nu }_2,{\nu }_1^0).$$

$${\dot{\mathbf{x}}}_i=\mathbf{F}({\mathbf{x}}_i)+K\sum _{j=1}^N{a}_{ij}[\mathbf{H}({\mathbf{x}}_j)-\mathbf{H}({\mathbf{x}}_i)],$$

$$R\dot{=}\frac{{\lambda }_N}{{\lambda }_2}<\frac{{\gamma }_u}{{\gamma }_l},$$

$${\lambda }_N(\alpha )={\nu }_{N_1}+\alpha {\mu }_{N_2}{\mathbf{v}}^{N_{1}T}{D}_1{\mathbf{v}}^{N_1}.$$

$${\lambda }_N(\alpha )=\alpha {\mu }_{N_2}+{\nu }_n^{\cancel{0}}.$$

$$L_{ϵ}({\mu }_i)=L_1+ϵ{\mu }_i{D}_1.$$

$${\lambda }_j(ϵ)$$

$${\mathbf{w}}^j(ϵ)$$

$${\mu }_i{D}_1{\mathbf{v}}^j+L_1{\hat{\mathbf{w}}}^j={\hat{\lambda }}_j{\mathbf{v}}^j+{\nu }_j{\hat{\mathbf{w}}}^j.$$

$${\hat{\lambda }}_j={\mu }_i{\mathbf{v}}^{jT}{D}_1{\mathbf{v}}^j$$

$${\lambda }_j(\alpha )={\nu }_j+\alpha {\mu }_i{\mathbf{v}}^{jT}{D}_1{\mathbf{v}}^j,$$

$${\alpha }^{-1}{L}_{ϵ}({\mu }_i)={\mu }_i{D}_1+ϵL_1,$$

$${\tilde{\lambda }}_j(ϵ)$$

$${\mathbf{w}}^j(ϵ)$$

$${\mu }_i{D}_1{\hat{\mathbf{w}}}^j+L_1\mathbf{x}={\hat{\lambda }}_j\mathbf{x}+{\mu }_i{\hat{\mathbf{w}}}^j.$$

$$\displaystyle\begin{array}[]{rl}\mu_{i}\hat{\bm{w}}^{\cancel{0}}+L_{1}^{\cancel{0}}\bm{x}^{\cancel{0}}&=\hat{\lambda}_{j}\bm{x}^{\cancel{0}}+\mu_{i}\hat{\bm{w}}^{\cancel{0}},\\ \to\hskip 8.61108ptL_{1}^{\cancel{0}}\bm{x}^{\cancel{0}}&=\hat{\lambda}_{j}\bm{x}^{\cancel{0}},\end{array}$$

$${\lambda }_j(\alpha )=\alpha {\mu }_i+{\nu }_j^{\cancel{0}}.$$

$${\tilde{\lambda }}_j(ϵ)$$

$${\mathbf{w}}^j(ϵ)$$

$${\mu }_i{D}_1{\hat{\mathbf{w}}}^j+L_1\mathbf{y}={\hat{\lambda }}_j\mathbf{y}.$$

$$L_1^0{\mathbf{y}}^0={\hat{\lambda }}_j{\mathbf{y}}^0,$$

Full text

I Introduction

The underlying structures that dictate the patterns of interactions that take place throughout nature and society are often described by complex networks Strogatz2001Nature . Examples of such networks include electrical power grids Motter2013NatPhys , faculty hiring networks Clauset2015SciAdv , gene regulatory networks Lee2002Science , and the structure of academic institutions Wang2017ANS . Many large networks are comprised of smaller subnetwork structures, for example motifs Milo2002Science , communities Girvan2002PNAS , layers DeDomenico2013PRX , self-similar structures Guimera2003PRE , or other subnetwork structures Gao2012NatPhys . In many such cases the properties of the larger network depends on properties of these smaller structures Chauhan2009PRE . Moreover, the collective organization of these smaller subnetwork structures often has a strong effect on the properties of many dynamical processes such as diffusion Gomez2013PRL , synchronization Skardal2012PRE and epidemic spreading Liu2005EPL .

Recently, Barrière et al. introduced the hierarchical product Barriere2009DAM1 ; Barriere2009DAM2 as a tool for building a large network using two smaller subnetworks. The hierarchical product is a generalization of the Cartesian product Hammack2011 , combining subnetworks in a less regular manner, resulting in a more disordered and heterogeneous structure – an important characteristic of many real-world networks Newman2003SIAM . Since its introduction, several properties of hierarchical products have been studied, including properties such as radius and diameter, clustering coefficient and degree distribution Barriere2016JPA . Recently, we provided an asymptotic analysis for the full spectrum of the adjacency matrix of the hierarchical product Skardal2016PRE . Here we apply these results in order to study the dynamics that take place on hierarchical products, particularly diffusion and synchronization.

Diffusion and synchronization represent two classical and well-studied classes of dynamical processes on networks. Diffusion has proven to be a particularly versatile tool in network science, identifying structural properties Lambiotte2011PRE ; Rosvall2014NatComm and serving as a mathematical model for other relaxation processes Grabow2012PRL ; Grabow2015PRE ; Skardal2016PREb . Synchronization dynamics are also strongly intertwined with the structures on which they evolve Pecora1998PRL , revealing topological properties of the the underlying networks Arenas2006PRL ; Pecora2014NatComm . The long-term dynamics of both diffusion and synchronization dynamics are determined by the eigenvalue spectrum of the network’s associated Laplacian matrix. In this work we apply techniques previously used to describe the eigenvalue spectrum of the adjacency matrix of hierarchical products to the Laplacian matrix in order to study the long-term diffusion and synchronization dynamics in this context. Our results allows us to extract the contributions that the two different subnetworks of the hierarchical product have on the long-term the macroscopic dynamics with respect to a coupling parameter that weighs of the secondary subnetwork relative to the the primary subnetwork. In the case of diffusion, two regimes emerge. For small coupling the diffusion rate is slow, increasing along with the coupling, and is determined by the structure of the secondary subnetwork. For large coupling the diffusion rate saturates and is determined by the structure of the primary subnetwork. Thus, a transition in both the dynamics and the contributions from the two subnetworks occurs as the coupling passes through this intermediate range. In the case of synchronization, three regimes emerge. For both small and large coupling the hierarchical product has poor synchronization properties, owing to a large deviation in a large gap between eigenvalues of the Laplacian that results in one subnetwork being weighted significantly more than the other. Synchronization properties are instead optimized at an intermediate, critical coupling value. Interestingly, this critical coupling value highlights the difference in overall connectivity in the primary and secondary subnetworks. Specifically, the critical coupling value is tuned to compensate for either the primary or secondary subnetwork having significantly weaker connectivity than the other. More broadly, the phenomena described in this paper identify the roles that primary and secondary subnetworks play in shaping large-scale dynamics in hierarchical products and which substructure promote vs hinder diffusion and synchronization.

The remainder of this paper is organized as follows. In Sec. II we define the hierarchical product and present asymptotic results describing the eigenvalue spectrum of the associated Laplacian matrix. In Sec. III we study the long-term behavior of diffusion on hierarchical products. We characterize the behavior of the smallest nontrivial eigenvalue, which dictates the diffusion rate and timescale. Using these results we identify two different regimes where diffusion processes behave differently. In Sec. IV we study the long-term behavior of synchronization on hierarchical products. Here we characterize both the largest and smallest nontrivial eigenvalues, which determine the synchronizability ratio. This allows us to identify regions of poor synchronization properties for small and large coupling, and optimal synchronization properties at an intermediate, critical coupling value. In Sec. V we conclude with a discussion of our results.

II The Hierarchical Product

II.1 Coupling Matrices

The hierarchical product represents a tool for building a large network from two smaller subnetworks. Here we will consider the hierarchical product of a primary and secondary subnetwork, denoted $G_{1}$ and $G_{2}$, respectively, each consisting of $N_{1}$ and $N_{2}$ nodes. We will assume that both networks are undirected and binary (or unweighted) so that the adjacency network $A_{1}$ and $A_{2}$ have entries $a_{ij}=a_{ji}=1$ if a link exists between nodes $i$ and $j$, and otherwise $a_{ij}=a_{ji}=0$. These properties can be generalized: in a directed network $a_{ij}$ and $a_{ji}$ need not be equal and in a weighted network $a_{ij}$ may take on values aside from zero and one. The hierarchical product of $G_{1}$ and $G_{2}$, denoted $G_{1}(U)\sqcap G_{2}$, is a network of $N=N_{1}\cdot N_{2}$ nodes that consists of $N_{2}$ copies of $G_{1}$ that are each connected to one another through the nodes in the root set $U$ via the topology of $G_{2}$. In Fig. 1 the hierarchical product is illustrated using an example with subnetworks $G_{1}$ and $G_{2}$ with $N_{1}=5$ and $N_{2}=4$ nodes, respectively. As a root set we use $U=\{1,4\}$, indicating that the four copies of $G_{1}$ are each connected by $G_{2}$ through nodes $1$ and $4$, so that in total $G_{1}(U)\sqcap G_{2}$ has $N=20$ nodes. The hierarchical product can be further generalized to include the product of an arbitrary number of subnetworks Barriere2009DAM2 , however such cases can be defined recursively, so we focus on hierarchical products of two subnetworks.

The goal of this work is to determine the effects that the two different subnetworks that make up a hierarchical product have on the long-term dynamics of diffusion and synchronization processes. To this end, we introduce a coupling parameter, denoted $\alpha$, that weighs the contribution of the secondary subnetwork $G_{2}$ in comparison to the primary subnetwork $G_{1}$. We incorporate this coupling parameter into the adjacency matrix of the hierarchical product, which is given by

[TABLE]

where $\otimes$ denotes the Kronecker product, $I_{2}$ is the $N_{2}\times N_{2}$ identity matrix, and $D_{1}$ is the $N_{1}\times N_{1}$ diagonal matrix whose $i^{\text{th}}$ diagonal entry is equal to one if vertex $i$ is in the root set $U$ and zero otherwise otherwise. Thus, $\alpha<1$ and $\alpha>1$ correspond to the links of the secondary subnetwork being weighted weaker and stronger, respectively, than the links of the primary subnetwork. While the spectral properties of Eq. (1) were studied in Ref. Skardal2016PRE , in this work we are interested in the Laplacian matrix due to the role it plays in dynamical processes, specifically diffusion and synchronization. For a network with adjacency matrix $A$, the Laplacian $L$ has entries defined $l_{ij}=\delta_{ij}\left(\sum_{j=1}^{N}a_{ij}\right)-a_{ij}=\delta_{ij}k_{i}-a_{ij}$, where $k_{i}$ denotes the nodal degree. In the case of a hierarchical product with adjacency matrix as in Eq. (1), the Laplacian matrix is given by

[TABLE]

where $L_{1}$ and $L_{2}$ are the Laplacian matrices of networks $G_{1}$ and $G_{2}$, respectively. The eigenvalue spectrum of the Laplacian matrix $L$ of any connected and undirected network has several important properties. First, since every row sums to zero there exists a trivial eigenvalue $\lambda_{1}=0$ whose associated eigenvector is constant, $\bm{w}^{1}\propto\bm{1}=[1,\dots,1]^{T}$. All other eigenvalues are real and positive, so they can be ordered $0=\lambda_{1}<\lambda_{2}\leq\dots\leq\lambda_{N}$. Finally, the eigenvectors are orthogonal and can therefore be normalized to form an orthonormal basis for $\mathbb{R}^{N}$.

II.2 Eigenvalues

The long-term dynamics of both diffusion and synchronization processes depend on the eigenvalues of the associated Laplacian matrix. Thus, for a full understanding of the long-term dynamics on the hierarchical product, we require a characterization of the eigenvalues associated with Eq. (2). In a previous publication we provided an asymptotic analysis for the eigenvalue spectrum of the adjacency matrix of the hierarchical product, i.e., Eq. (1) Skardal2016PRE . Here we apply the same methodology to the case of the Laplacian. Our goal is to classify the eigenvalues of $L_{\alpha}$ in terms of the eigenvalues and eigenvectors of the Laplacian of its subnetworks $L_{1}$ and $L_{2}$, as well as the coupling parameter $\alpha$ and the root set encapsulated in the matrix $D_{1}$. We will denote the eigenvalues of $L_{1}$ and $L_{2}$ as $\nu_{i}$ and $\mu_{i}$, respectively, and the associated eigenvectors $\bm{v}^{i}$ and $\bm{u}^{i}$, respectively. We seek the eigenvalues, denote $\lambda$, of $L_{\alpha}$. Before classifying them, we illustrate the general behavior of the eigenvalues of the hierarchical product as a function of the coupling parameter $\alpha$ in Fig 2, using the network illustrated in Fig. 1 as an example. Broadly speaking, these eigenvalues split into two groups for both small and large $\alpha$. For small $\alpha$ one group of eigenvalues are themselves small, scaling approximately as $\alpha$ and the other remains approximately constant, on the order of one. For large $\alpha$ one group of eigenvalues also remains approximately constant, on the order of one, but the other group is itself large, also scaling approximately as $\alpha$. When $\alpha$ is itself on the order of one these groups coalesce in a complicated arrangement.

The classification of the eigenvalues $\lambda$ then begins with the analysis of a new set of matrices. Specifically, $\lambda$ is an eigenvalue of $L_{\alpha}$ if and only if it is also an eigenvalue of one of the matrices given by

[TABLE]

where $L_{\alpha}(\mu_{i})$ is one of the $N_{2}$ possible $N_{1}\times N_{1}$ matrix constructed via a linear combination of $L_{1}$ and $D_{1}$, where $D_{1}$ is scaled by one of the $N_{2}$ eigenvalues $\mu_{i}$ of $L_{2}$ Barriere2009DAM1 . In total, there are $N_{2}$ such matrices $L_{\alpha}(\mu_{i})$, each of which have $N_{1}$ eigenvalues, resulting in the full spectrum of $N_{1}\cdot N_{2}$ eigenvalues of $L_{\alpha}$. Thus, the eigenvalue problem of $L_{\alpha}$ is reduced to the set of smaller eigenvalue problems for the collection of $L_{\alpha}(\mu_{i})$.

While the collection of eigenvalues of Eq. (3) can be found perturbatively as in Ref. Skardal2016PRE , a specific subset deserves particular attention here. Recall that the Laplacian of a connected network has precisely one zero eigenvalue. Thus, $L_{2}$ has one zero eigenvalue $\mu_{1}=0$. Inserting this into Eq. (3) yields, simply,

[TABLE]

Eq. (4) implies then that precisely $N_{1}$ of the eigenvalues of $L_{\alpha}$ are given by the eigenvalues $\nu_{i}$ of $L_{1}$, and that these eigenvalues remain constant regardless of the value of $\alpha$. Exactly one of these eigenvalues corresponds to the zero eigenvalue $\nu_{1}=0$ of $L_{1}$, with the other $N_{1}-1$ being finite.

As for the remaining $N_{1}(N_{2}-1)$ eigenvalues of $L_{\alpha}$, we apply a perturbative analysis for the limits of small and large coupling. Here we present the results and leave the details for the interested reader in Appendix A. In the case of small coupling, i.e., $\alpha\ll 1$ the eigenvalues are given, to first order in $\alpha$, by

[TABLE]

which represents the contributions of the two subnetworks via $\nu_{j}$ and $\mu_{i}$ ($j=1,\dots,N_{1}$, $i=2,\dots,N_{2}$) to the eigenvalue spectrum. It should be noted that using $\mu_{1}=0$ in Eq. (5) recovers the constant eigenvalues from Eq. (4), however we consider this a separate case since the eigenvalues from Eq. (4) are exact, while those in Eq. (5) are approximations.

In the limit of large coupling, i.e., $\alpha\gg 1$, the analysis becomes more complicated with the $N_{1}(N_{2}-1)$ eigenvalues splitting into two distinct groups due to the degeneracy of $D_{1}$. The first group yields $n(N_{2}-1)$ eigenvalues, where $n$ is the size of the root set and the number of nonzero entries of $D_{1}$, and are given by

[TABLE]

where $\nu_{j}^{\cancel{0}}$ is an eigenvalues of the $n\times n$ matrix $L_{1}^{\cancel{0}}$ constructed by removing all rows and columns of $L_{1}$ corresponding to zero diagonal entries of $D_{1}$. The remaining $(N_{1}-n)(N_{2}-1)$ eigenvalues are then given by

[TABLE]

where $\nu_{j}^{0}$ is an eigenvalues of the $(N_{1}-n)\times(N_{1}-n)$ matrix $L_{1}^{0}$ constructed by removing all rows and columns of $L_{1}$ corresponding to nonzero diagonal entries of $D_{1}$. In Figs. 3(a) and (b) we compare the analytical approximations (dashed blue and dot-dashed red curves) for the eigenvalues of $L_{\alpha}$ to the actual values (solid black) for the example hierarchical product illustrated in Fig. 1 for small and large $\alpha$. In both cases the approximations accurately describe the behavior of the eigenvalues for sufficiently small and large $\alpha$, with the approximations breaking down when $\alpha$ is approximately of order one.

III Diffusion

We now turn our attention to the long-term dynamics of diffusion process on hierarchical products. Given an adjacency matrix with entries $a_{ij}$, diffusion is governed by the following equations

[TABLE]

which can be rewritten in vector form as

[TABLE]

where $\bm{x}$ is the state vector of the process and $L$ is the Laplacian. We note here that we focus on the specific case of diffusion related to heat transfer and relaxation dynamics, in which case we use the combinatorial Laplacian $L=D-A$, where $D=\text{diag}(n_{1}\dots,k_{N})$. In the case of diffusion related to a random walk processes, the symmetric or asymmetric versions of the normalized Laplacian, $L=I-D^{-1/2}AD^{-1/2}$ or $L=I-D^{-1}A$, may be used. We note that in either case the methodology for approximating eigenvalues may be preserved, but in the asymmetric case the emergence of complex eigenvalues may require more care when applying these results.

Assuming that the underlying network is connected, the dynamics of Eqs. (8) and (9) relax to the steady state $x_{1}=\cdots=x_{N}=x_{\infty}$ in the limit $t\to\infty$. This relaxation is exponential, specifically with

[TABLE]

i.e., the rate of diffusion is given by smallest nontrivial eigenvalue $\lambda_{2}$ and the timescale of diffusion is given by its inverse $\lambda_{2}^{-1}$. Therefore, we seek specifically the smallest nontrivial eigenvalue $\lambda_{2}$ from our approximations above. In the small coupling regime, $\alpha\ll 1$, the full set of eigenvalues is given by the collection of $N_{1}$ eigenvalues of $L_{1}$ along with the $N_{1}(N_{2}-1)$ eigenvalues in Eq. (5). Since the nontrivial eigenvalues of $L_{1}$ are all of order one, the smallest nontrivial eigenvalue is given by using $\nu_{1}=0$ and $\mu_{2}$ in Eq. (5), resulting in

[TABLE]

where we have used that $\bm{v}_{1}=\bm{1}/\sqrt{N_{1}}$. In the large coupling regime, $\alpha\gg 1$ the full set of eigenvalues is given, again, by the collection of $N_{1}$ eigenvalues of $L_{1}$, along with those given in Eqs. (6) and (7). Inspecting all possible combinations, the smallest nontrivial eigenvalue is then given by either the smallest nontrivial eigevnalue of $L_{1}$, or the smallest eigenvalue of $L^{0}$, i.e.,

[TABLE]

In general, it is impossible to determine a priori which eigenvalue in Eq. (14) in smallest; as we shall see below, different combinations of different subnetwork structures yield different outcomes.

We first compare our the predictions of our approximations to direct simulation results. In Fig. 4 we plot our theoretical prediction of the diffusion rate (dashed black), i.e., Eqs. (13) and (14) for the small and large coupling regimes, respectively, to the diffusion rate observed from simulations (blue circles) on the hierarchical product illustrated in Fig. 1. Simulated results are computed by fitting the simulations after a significant transient to an exponential. The diffusion timescale is plotted in the inset. The two different dynamical behaviors, corresponding to small and large coupling, are observed and are well captured by the predictions. Specifically, in the small coupling regime the diffusion rate is very slow, and scales with the coupling parameter $\alpha$, which can be observed directly from Eq. (13). Moreover, Eq. (13) reveals that the long-time diffusion dynamics in the small coupling regime is completely determined by the structure of the secondary subnetwork via the eigenvalue $\mu_{2}$ and the size of the root set via the fraction $n/N_{1}$.

In contrast to the small coupling regime, in the large coupling regime, the diffusion rate saturates to the order one value given in Eq. (14). Moreover, Eq. (14) reveals that in this regime the long-time diffusion dynamics in the large coupling regime are completely determined by the structure of the primary subnetwork via the eigenvalue $\nu_{2}$, and possibly in combination with the root set via the eigenvalue $\nu_{1}^{0}$. This highlights a transition between small and large coupling from two perspectives. First, the rate of diffusion itself is increasing, scaling with $\alpha$, for small coupling, and saturates to a constant value for large $\alpha$. Second, the role that the components of the hierarchical product play in this behavior changes; for small $\alpha$ dynamics are dictated by the secondary subnetwork and for large $\alpha$ dynamics are dictated by the primary subnetwork.

Next we investigate in more detail the large coupling regime specifically the determination of $\lambda_{2}$ as $\nu_{2}$ or $\nu_{1}^{0}$ in Eq. (14). In both cases we note that the structure of the secondary subnetwork is irrelevant – only the primary subnetwork and the root set determine these quantities. However, in which cases $\lambda_{2}$ is determined by with choice in unclear. To better understand these quantities, we compare them in Fig. 5, plotting $\nu_{2}$ (blue circles) and $\nu_{1}^{0}$ (red triangles) computed from a collection of $1000$ realizations of Erdős-Rényi (ER) networks Erdos1960 of size $N=100$ constructed with link probability $p=0.1$ as a function of different root sets fractions $n/N_{1}$, where nodes in the root sets are randomly chosen. Results represent the average over the $1000$ networks, with the standard deviation denoted by dashed curves. As the root set fraction increases the value $\nu_{2}$ remains constant (as should be expected for a set network model) and $\nu_{1}^{0}$ increases from zero. Thus, for small $n/N_{1}$ we have $\nu_{1}^{0}<\nu_{2}$, indicating that $\lambda_{2}=\nu_{1}^{0}$, but for large enough $n/N_{1}$ we have $\nu_{2}<\nu_{1}^{0}$, indicating that $\lambda_{2}=\nu_{2}$. For the network model chosen here we find that this transition occurs at $n/N_{1}\approx 0.3487$, which is illustrated with the vertical dot-dashed line. Physically, this suggests that for a small enough root set the diffusion rate is determined by a combination of the structure of the primary subnetwork and the root set itself, but for a large enough root set the diffusion rate is determined solely by the primary subnetwork.

Given that the size of the root set plays a role in whether $\lambda_{2}$ is given by $\nu_{2}$ or $\nu_{1}^{0}$, it is natural to ask whether the particular locations of the root set also plays a role. In other words, does the behavior of $\nu_{2}$ and $\nu_{1}^{0}$ depend significantly on which nodes belong to the root set, in addition to the size of it? The ER model used in Fig. 5 yields relatively homogeneous networks where nodes have by-and-large very similar structural properties. To investigate this new question we then use the Barabasi Albert (BA) model Barabasi1999Science , which yields much more heterogeneous networks. Specifically, we consider BA networks of size $N=100$ with minimum degree $k_{0}=3$, but choose the root set in three different ways. In Fig. 6 we plot the average behavior of $\nu_{2}$ (blue circles) and $\nu_{1}^{0}$ (red triangles) as a function of the root set fraction $n/N_{1}$, choosing the root set to contain (a) the highest degree nodes in the network, (b) randomly selected nodes, and (c) the lowest degree nodes in the network. The generic behavior is similar in the sense that $\nu_{2}$ remains constant and $\nu_{1}^{0}$ increases with $n/N_{1}$. However, the critical root set fraction at which the transition from $\lambda_{2}=\nu_{1}^{0}$ to $\lambda_{2}=\nu_{2}$ occurs is different. Specifically, when the root set consists of the highest degree nodes in the network this transition occurs quite early, at $n/N_{1}\approx 0.1630$. Conversely, when the root set consists of the lowest degree nodes in the network this transition occurs quite late, at $n/N_{1}\approx 0.3843$. When the root set consists of randomly chosen nodes this transition occurs in between, at $n/N_{1}\approx 0.3375$. Thus, when the root set consists of lower degree nodes, it plays a role in determining the diffusion rate for a larger range of root set fractions than when it consists of higher degree nodes.

IV Synchronization

Next we turn to synchronization on hierarchical products. Specifically, we consider synchronization of identical, chaotic dynamical systems, whose dynamics are governed by

[TABLE]

where $\bm{x}_{i}$ is the state vector of node $i$, $\bm{F}(\bm{x})$ is the (assumed chaotic) vector field describing the internal dynamics of each node, $K$ is the global coupling strength, and $\bm{H}(\bm{x})$ is the coupling function. The dynamics of Eq. (15) are typically treated by studying the stability of the synchronized state $\bm{x}_{1}(t)=\cdots=\bm{x}_{N}(t)$, which can be determined using the Master Stability Function (MSF) approach Pecora1998PRL . In particular, the synchronized state is linear stable if all the nontrivial eigenvalues of the Laplacian matrix scaled by the coupling strength $K$ fall within an appropriately defined region of stability. (For the sake of brevity, we forgo a discussion of further technical details and refer the interested reader to the original work in Ref. Pecora1998PRL .) While the region of stability depends on the particular dynamical system $\bm{F}$ and the coupling function $\bm{H}$ in Eq. (15), in many cases the region of stability is a finite interval, denoted $[\gamma_{l},\gamma_{u}]$ Huang2009PRE . Synchronization can then be achieved if a coupling $K$ can be chosen such that $K\lambda_{i}$ for $i=2,\dots,N$ fall within the interval. This is true if and only if the eigenvalues satisfy

[TABLE]

where $R$ is the synchronizability ratio of the network. In particular, the smaller the synchronizability ratio $R$ a given network has, the more synchronizable it is.

The synchronizability of a given hierarchical product then requires both the smallest and largest nontrivial eigenvalues, $\lambda_{2}$ and $\lambda_{N}$. Since we characterized in detail the smallest nontrivial eigenvalue in the previous section, we now turn to the largest. In the small coupling regime we refer back to Eq. (5). Note that these eigenvalue are always larger than the corresponding constant eigenvalues of $L_{1}$ since $\alpha\mu_{i}\bm{v}^{jT}D_{1}\bm{v}^{j}>0$. Thus, the maximum is obtained by choosing $j=N_{1}$ and $i=N_{2}$, yielding

[TABLE]

In the large coupling regime, we find that the largest eigenvalue is given by Eq. (6), using $i=N_{2}$ and $j=n$, yielding

[TABLE]

In Figs. 7(a) and (b) we demonstrate how the synchronizability of the hierarchical product illustrated in Fig. 1 behaves as a function of the coupling parameters, first plotting the separate behaviors of the actual (solid black) approximate (dashed blue and dot-dashed red) values of $\lambda_{N}$ and $\lambda_{2}$ in panel (a), then in panel (b) the actual (solid black) and approximate (dashed blue) synchronizability ratio $R=\lambda_{n}/\lambda_{2}$. Specifically, we see in panel (b) that for both very large and very small $\alpha$ the synchronizability ratio is large, indicating that the hierarchical product has poor synchronization properties. This is due to the large gap between $\lambda_{2}$ and $\lambda_{N}$ which can be observed directly in panel (a), and can be physically attributed to one of the two subnetworks being weighted much heavier than the other. Instead, the hierarchical product displays the best synchronizability ratio for intermediate values of $\alpha$, suggesting that hierarchical products have have the best synchronization properties when the two subnetworks are weighted roughly equally.

Next we investigate the role that the two different subnetworks and the root set play in determining the sychronizability of the hierarchical product. First we consider the synchronizability itself – specifically the optimal (minimal) synchronizability attainable for a given hierarchical product. We find that this quantity depends primarily on the size of the root set. In Fig. 8 we plot the actual (blue circles) and approximate (red triangles) optimal synchronizability ratio $R_{\text{min}}$ vs. the root set fraction $n/N_{1}$ found for hierarchical products constructed using ER networks of size $N_{1}=50$ and $N_{2}=20$ using link probabilities $p=0.2$ and $0.5$, respectively. (These probabilities are chosen to attain a rough balance of the mean degree.) Results represent an average over $100$ networks, with standard deviation denoted by dashed curves. In general we see that the larger the root set fraction, the more synchronizable the hierarchical product can be when $\alpha$ is properly tuned. Thus, the more pathways built into the hierarchical product via the root set, the more favorable the synchronization properties.

A more interesting question, however, is at what coupling value $\alpha_{c}$ is the synchronizability of a hierarchical product optimized? We find that this critical coupling value does not depend significantly on the root set itself, but rather the contrast between the primary and secondary subnetworks. In fact, the answer to this question sheds light on role of the primary and secondary networks in relation to one another. We consider hierarchical products constructed from ER networks, both of size $N_{1}=N_{2}=50$, for each subnetwork choosing the link probabilities randomly to allow the mean degrees for the subnetworks, denoted $\langle k\rangle_{1}$ and $\langle k\rangle_{2}$, to vary between $5$ and $45$. Using the mean degree of each subnetwork as a proxy for overall connectivity, we then compare the critical coupling values $\alpha_{c}$ to the connectivity ratio $\langle k\rangle_{1}/\langle k\rangle_{2}$, plotting in Fig. 9 the results from $100$ different networks the actual (blue circles) and approximated (red triangles) results. Figure 9 shows a positive, roughly power-law relationship between $\alpha_{c}$ and $\langle k\rangle_{1}/\langle k\rangle_{2}$. This suggests that, to achieve optimal synchronizability, the coupling should be tuned to balance the connectivity properties of the primary and secondary subnetworks. If the ratio $\langle k\rangle_{1}/\langle k\rangle_{2}$ is large (i.e., larger than one), indicating that the primary subnetwork is more strongly connected than the secondary subnetwork, then the coupling should be increased to strengthen the secondary subnetwork in compensation. On the other hand, If the ratio $\langle k\rangle_{1}/\langle k\rangle_{2}$ is small (i.e., less than one), indicating that the primary subnetwork is connected more weakly than the secondary subnetworks, then the coupling should be decreased to weaken the secondary subnetwork in compensation. Moreover, we find that the results fall roughly around the power-law relationship $\alpha_{c}\propto\left(\langle k\rangle_{1}/\langle k\rangle_{2}\right)^{\beta}$ for $\beta\approx 1.18$, as illustrated with the dashed black curve.

V Discussion

In this paper we have studied the long-term dynamics of diffusion and synchronization processes on hierarchical products. We have applied the methodology from previous work Skardal2016PRE characterizing the eigenvalues of the adjacency matrix of hierarchical products to the eigenvalues of the Laplacian matrix, allowing us to make analytical predictions for both diffusion and synchronization dynamics. In particular, this has allowed us to identify the roles that the primary and secondary subnetworks play in shaping the long-term dynamics in relation to a coupling parameter that weighs the contribution of the secondary subnetwork relative to the primary subnetwork. More generally our results explore the effects that different substructure of networks play in shaping large-scale dynamics by either promoting or inhibiting these processes.

In the case of diffusion, we have identified two regimes corresponding to small and large coupling. In the small coupling regime the diffusion rate is slow, scaling with the coupling itself, and is completely determined by the structure of the secondary subnetwork. In the large coupling regime the diffusion rate saturates to a constant value which is determined by the structure of the primary subnetwork, as well as possibly the size of the root set and its structure. Thus, there is an transition that occurs as coupling is varied through intermediate values, both in terms of the long-term dynamical behavior, as well as the roles that the different structures that make up the hierarchical product play in shaping those dynamics.

In the case of synchronization, we find that the synchronization properties of hierarchical products are poor in both the small and large coupling regimes, but is optimized at an intermediate critical coupling value that minimizes the synchronizability ratio. In general, the optimal synchronizability ratio that a hierarchical product can attain, assuming $\alpha$ can be properly tuned, improves as the size of the root set increases. However, a more interesting phenomenon occurs with the critical coupling parameter that optimizes synchronization, which highlights the difference in overall connectivity between the primary and secondary subnetworks. Specifically, the critical coupling value is tuned to compensate for this difference, either strengthening or weakening the secondary subnetwork to bring is connectivity closer to that of the primary subnetwork.

Throughout this work we have focused on the case of undirected subnetworks, resulting in undirected hierarchical product. In the case of directed subnetworks it is straightforward to see that the resulting hierarchical product also becomes directed. In principle, the techniques used here to calculate eigenvalues may be preserved, however the emergence of complex eigenvalues may require some care when applying these results. This is also true when working with the asymmetric normalized Laplacian matrix for random walks, even in the case of undirected networks.

The class of networks investigated here, i.e., the hierarchical product Barriere2009DAM1 ; Barriere2009DAM2 , represents a relatively wide subset possible generalizations of classical graph products Hammack2011 . In general, graph products represent natural ways of building larger networks from two or more smaller subnetworks where the macroscopic properties of the large network can be understood in terms of the properties of the smaller subnetworks that comprise it. The general notion of a network consisting of smaller substructures remains a central theme in physics and mathematics, with examples including multilayer and multiplex networks Granell2013PRL ; Kouvaris2015SciRep ; DeDomenico2016NatPhys ; Nicosia2017PRL , modular networks Liu2005EPL ; Skardal2012PRE , hierarchical and hierarchical modular networks Boettcher2008JPA ; Boettcher2009PRE ; Kaiser2010Front ; Moretti2013NatComms , and networks of networks Gao2011PRL ; Bianconi2014PRE . Many of these cases share commonalities, for example the behavior of the Laplacian eigenvalues we observe in hierarchical products (e.g., see Figs. 2 and 4) reflects the behavior of Laplacian eigenvalues in multiplex networks Gomez2013PRL ; Sole2013PRE . Given this overlap in phenomenological behavior, we hypothesize that understanding the macroscopic structural properties of hierarchical products is not only important in the context of graph products, but also more generally for wider classes of networks. To date, a handful of studies have investigated the overall structure of hierarchical products Barriere2016JPA ; Skardal2016PRE , little work has focused on behavior of dynamical processes taking place on hierarchical products.

Finally, we emphasize that the contributions of this paper, i.e., the description of the long-term diffusion and synchronization dynamics on hierarchical products, fit within the broader question of how various structures and organizations in complex networks dictate large-scale dynamical processes. Specifically, the findings presented here can be interpreted as investigating how different components and substructures of a given network, and their relative strengths, function in shaping the dynamics that occur across the whole network. This broad question has been investigated for various kinds of networks (e.g., see those listed above); here we study this broad question in the context of a graph product. In particular, the long-term behaviors of both diffusion and synchronization dynamics identify the role of the secondary subnetwork as a connector in comparison to the more central primary subnetwork, as well as the role that nodes in the root set play in facilitating these connections. Therefore, the coupling parameter modifies the relative strengths of the overall connectivity within the primary subnetworks compared to the connectivity between different subnetworks. Specifically, it is the secondary subnetwork structure that is responsible for the smaller eigenvalues in the small coupling regime, whereas in the large coupling regime the primary subnetwork is responsible for the smaller eigenvalues. More broadly, the transitions that we observe in the dynamics can be interpreted as a shift in which of these connectivities becomes the effective bottleneck for the dynamics and the primary hinderance for diffusion or spread of consensus (i.e., synchronization) throughout the network.

Appendix A Eigenvalue Perturbation Analysis

Here we present the perturbative analysis for the eigenvalues of the Laplacian $L_{\alpha}$ in Eq. (2), which are in turn given by the eigenvalues of Eq. (3). We consider here the $N_{1}(N_{2}-1)$ eigenvalues corresponding to inserting the nonzero eigenvalues $\mu_{2},\dots,\mu_{N_{2}}$ into Eq. (3). Beginning with the limit of small coupling, $\alpha\ll 1$, we proceeding perturbatively as in Ref. Skardal2016PRE , we make the common notational change $\epsilon=\alpha$ such that $\epsilon\ll 1$ is a small parameter and study the eigenvalues of

[TABLE]

In the limit $\epsilon\to 0^{+}$ we recover the spectrum of $L_{1}$, i.e., eigenvalues $\nu_{i}$ and eigenvectors $\bm{v}^{i}$, and therefore we propose a perturbative ansatz of the form

[TABLE]

and seek the coefficient $\hat{\lambda}_{j}$ of the first-order correction, i.e., searching for the leading order behavior of the Taylor series for $\lambda_{j}(\epsilon)$ and $\bm{w}^{j}(\epsilon)$. Inserting Eqs. (19), (20), and (21) into the eigenvalue equation $L_{\epsilon}(\mu_{i})\bm{w}^{j}(\epsilon)=\lambda_{j}(\epsilon)\bm{w}^{j}(\epsilon)$ and collecting the leading order terms at $\mathcal{O}(\epsilon)$, we obtain

[TABLE]

Left-multiplying Eq. (22) by $\bm{v}^{jT}$ and noting that the term on the left-hand side $\bm{v}^{j}L_{1}\hat{\bm{w}}^{j}=\nu_{j}\bm{v}^{j}\hat{\bm{w}}^{j}$ cancels with the right-hand side, we obtain

[TABLE]

We note that terms similar to the right hand side in Eq. (23) appear often in perturbaitve analyses, and are akin to the first order correction to the energy of a Hamiltonian Schrodinger . Substituting back $\epsilon=\alpha$, we have that the eigenvalues of $L_{\alpha}(\mu_{i})$ to leading order are given by

[TABLE]

giving the expression presented in Eq. (5) in the main text.

Next we consider the limit of large coupling, $\alpha\gg 1$, now letting $\epsilon=\alpha^{-1}$ be a small parameter. We again proceed perturbatively, noting that now

[TABLE]

so that after finding the eigenvalues of the right hand side of Eq. (25) for $\mu_{2},\dots,\mu_{N_{2}}$, we then multiply by $\alpha$ to obtain the final eigenvalues. As in Ref. Skardal2016PRE , the perturbative analysis for large coupling then becomes more complicated than that for small coupling due to the fact that when $\epsilon=0$ the right hand side of Eq. (25) reduces to the matrix $\mu_{i}D_{1}$, which is degenerate. Specifically, $D_{1}$ has precisely $n$ eigenvalues equal to one and $(N_{1}-n)$ eigenvalues equal to zero, where $n$ is the number of nonzero entries of $D_{1}$, i.e., the size of the roots set $U$. We will refer to the eigenspaces associated with the one and zero eigenvalues of $D_{1}$ as the nontrivial and trivial eigenspaces. (Note that the trivial eigenspace is precisely the nullspace.) Specifically, the nontrivial eigenspace of $D_{1}$ is the span of all vectors whose entries are zero where the diagonal entries of $D_{1}$ are zero, and the trivial eigenspace of $D_{1}$ is the span of all vectors whose entries are zero where the diagonal entries of $D_{1}$ are non-zero. This requires us to consider two subcases of our asymptotic analysis: one for the nontrivial eigenspace which will yield $n$ eigenvalues and another for the trivial eigenspace which will yield $N_{1}-n$ eigenvalues.

We begin with the nontrivial eigenspace of $D_{1}$ and propose a perturbative ansatz of the form

[TABLE]

where the vector $\bm{x}$ is in the non-trivial nullspace of $D_{1}$, i.e., $D_{1}\bm{x}=\bm{x}$. Inserting Eqs. (25), (26), and (27) into the eigenvalue equation $\alpha^{-1}L_{\epsilon}(\mu_{i})\bm{w}^{j}(\epsilon)=\tilde{\lambda}_{j}(\epsilon)\bm{w}^{j}(\epsilon)$ and collecting the leading order terms at $\mathcal{O}(\epsilon)$, we obtain

[TABLE]

Next, the entries of $\bm{x}$ that correspond to zeros in the diagonal of $D_{1}$ (i.e., nodes that do not belong to the root set $U$) are zero, so we eliminate these $(N_{1}-n)$ entries from Eq. (28) and obtain the following $n$-dimensional vector equation:

[TABLE]

where $L_{1}^{\cancel{0}}$ is the $n\times n$ matrix obtained by keeping the rows and columns of $L_{1}$ corresponding to non-zero diagonal entries of $D_{1}$ and similarly $\hat{\bm{w}}^{j\cancel{0}}$ and $\bm{x}^{\cancel{0}}$ are the $n$-dimensional vectors obtained by keeping the same entries of $\hat{\bm{w}}^{j}$ and $\bm{x}$. Thus, $\hat{\lambda}_{j}$ is one of the $n$ eigenvalues of the matrix $L_{1}^{\cancel{0}}$, which we will denote $\nu_{j}^{\cancel{0}}$. Inserting this back into Eq. (26), replacing $\epsilon=\alpha^{-1}$, and multiplying by $\alpha$, the $n$ eigenvalues of $L_{\alpha}(\mu_{i})$ corresponding to the nontrivial eigenspace of $D_{1}$ to leading order are given by

[TABLE]

Turning our attention to the trivial eigenspace of $D_{1}$, we introduce a new perturbative anstaz:

[TABLE]

where the vector $\bm{y}$ is now in the nullspace of $D_{1}$, i.e., $D_{1}\bm{y}=\bm{0}$. Inserting Eqs. (25), (33), and (34) into the eigenvalue equation $\alpha^{-1}L_{\epsilon}(\mu_{i})\bm{w}^{j}(\epsilon)=\tilde{\lambda}_{j}(\epsilon)\bm{w}^{j}(\epsilon)$ and collecting the leading order terms at $\mathcal{O}(\epsilon)$, we obtain

[TABLE]

Similarly to Eq. (28), several vector entries in Eq. (35) are zero: this time all entries of $\bm{y}$ that correspond to ones in the diagonal of $D_{1}$ (i.e., nodes that are in the root set $U$) are zero. We therefore eliminate these $n$ entries from Eq. (35) to obtain the following $(N_{1}-n)$-dimensional vector equation

[TABLE]

where $L_{1}^{0}$ is the $(N_{1}-n)\times(N_{1}-n)$ matrix obtained by keeping the rows and columns of $L_{1}$ corresponding to zero diagonal entries of $D_{1}$ and similarly $\bm{y}^{0}$ is the $(N_{1}-n)$-dimensional vector obtained by keeping the same entries of $\bm{y}$. Thus, $\hat{\lambda}_{j}$ is an eigenvalue of the matrix $L_{1}^{0}$, which we will denote $\nu_{j}^{0}$. Inserting this back into Eq. (33), replacing $\epsilon=\alpha^{-1}$, and multiplying by $\alpha$, the $(N_{1}-n)$ eigenvalues of $L_{\alpha}(\mu_{i})$ corresponding to the trivial eigenspace of $D_{1}$ to leading order are given by

[TABLE]

Eqs. (32) and (37) give those expressions presented in Eqs. (6) and (7) in the main text.