.

# Geometric unfolding of synchronization dynamics on networks

Lluís Arola-Fernández  Departament d’Enginyeria Informàtica i Matemàtiques,
Universitat Rovira i Virgili, 43007 Tarragona, Spain  Per Sebastian Skardal
Department of Mathematics, Trinity College, Hartford, CT 06106, USA  Alex
Arenas<EMAIL_ADDRESS>Departament d’Enginyeria Informàtica i
Matemàtiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain

###### Abstract

We study the synchronized state in a population of network-coupled,
heterogeneous oscillators. In particular, we show that the steady-state
solution of the linearized dynamics may be written as a geometric series whose
subsequent terms represent different spatial scales of the network. Namely,
each addition term incorporates contributions from wider network
neighborhoods. We prove that this geometric expansion converges for arbitrary
frequency distributions and for both undirected and directed networks provided
that the adjacency matrix is primitive. We also show that the error in the
truncated series grows geometrically with the second largest eigenvalue of the
normalized adjacency matrix, analogously to the rate of convergence to the
stationary distribution of a random walk. Lastly, we derive a local
approximation for the synchronized state by truncating the spatial series, at
the first neighborhood term, to illustrate the practical advantages of our
approach.

††preprint: AIP/123-QED

> Synchronization in ensembles of network-coupled heterogeneous oscillators
> plays a vital role in a number of natural and engineered phenomena, ranging
> from cell cycles to robust power systems Prindle  _et al._ (2012); Witthaut
> and Timme (2012). In a wide range of practical applications of such systems,
> including synchronization optimization Skardal, Taylor, and Sun (2014,
> 2016); Taylor, Skardal, and Sun (2016), control Skardal and Arenas (2015,
> 2016), and global entrainment Dörfler, Chertkov, and Bullo (2013), the
> particular structure of the the fully synchronized state or a target
> synchronized state is critical for achieving the task at hand. In this
> paper, we study the synchronized state in networks of heterogeneous
> oscillators and show that it can written as a particular geometric series
> where subsequent series terms describe contributions from larger network
> neighborhoods. This spatial description of the state provides a contrast to,
> e.g., spectral methods that rely on global information. In particular,
> truncation of the geometric series at a fixed order is equivalent to
> building the synchronized state using only local, i.e., decentralized,
> information. We provide a rigorous proof of our main result, show that error
> in the truncated series decays geometrically with the second largest
> eigenvalue of the normalized adjacency matrix, and illustrate the practical
> advantages of this approach.

##  I Introduction

The study of synchronization dynamics on networks of heterogeneous oscillators
has shed light on a wide range of natural and engineered phenomena and has
provided significant insights into structural properties of complex
networksWatts and Strogatz (1998); Arenas  _et al._ (2008); Pikovsky,
Rosenblum, and Kurths (2003). A particularly useful model for these purposes
uses phase oscillators that evolve according to the following equations of
motion,

$\dot{\theta}_{i}=\omega_{i}+\sigma\sum_{j=1}^{N}a_{ij}H(\theta_{j}-\theta_{i}),$
(1)

where $\theta_{i}$ is the phase of oscillator $i$, $i=1,\dots N$, $\omega_{i}$
is its natural frequency, $\sigma$ is the global coupling strength, and
$H(\theta_{j}-\theta_{i})$ is a $2\pi$-periodic non-linear coupling function
of the state differences of pairs of connected oscillators. The network
connectivity and the relative intensity of the interactions are captured by
the entries $a_{ij}$ of the $N\times N$, real-valued adjacency matrix $A$ that
in general can be undirected or directed, binary or weighted. While the
particular choice of $H(\theta)$ depends on the application at hand, notable
choices $H(\theta)=\sin(\theta)$ and $H(\theta)=\sin(\theta-\alpha)$, where
$\alpha\in[-\pi/2,\pi/2]$, yield the network Kuramoto and Sakaguchi-Kuramoto
models, respectively Kuramoto (2003a).

A particularly suitable framework useful for analyzing the interplay between
structure and dynamics, described here by the network $A$ and the natural
frequencies $\bm{\omega}$, respectively, emerges when a strongly synchronized,
phase-locked state becomes attainable, for which we assume that $|H(0)|\ll 1$,
so that for a strongly synchronized state where $|\theta_{j}-\theta_{i}|$ we
have that $H(\theta_{j}-\theta_{i})\approx
H(0)+H^{\prime}(0)(\theta_{j}-\theta_{i})$Skardal, Taylor, and Sun (2014). (We
note here that for the classical Kuramoto model we have $H(0)=0$ and
$H^{\prime}(0)=1$.) By defining the effective frequency
$\hat{\omega}_{i}=\omega_{i}+\sigma H(0)k_{i}^{in}$ and entering a suitable
rotating frame, the synchronized state is then described by a the fixed point
satisfying $\dot{\bm{\theta}}=0$. After linearizing, the synchronized state
then satisfies the following equation,

$\bm{\hat{\omega}}=\sigma H^{\prime}(0)L\bm{\theta},$ (2)

where $L=D-A$ is the Laplacian of the network, and $D$ is a diagonal matrix
with the in-degrees or in-strengths of the nodes, i.e.
$D_{ii}=k_{i}=\sum_{j=1}^{N}a_{ij}$. Since the matrix ${L}$ has zero row sum,
it has a trivial eigenvalue $\lambda_{1}=0$ with a constant associated
eigenvector $\bm{v}_{1}\propto\bm{1}$, and therefore it is singular and not
invertibleNewman (2010); Skardal, Taylor, and Sun (2014). Moreover, this
spectral property reveals an important physical characteristic of the system,
namely that the dynamics are invariant to a constant shift to the phases,
i.e., translation along the synchronization manifold defined as the span of
the trivial eigenvector $\bm{v}_{1}$. Thus, while solutions to Eq. (2) are not
unique, that which minimizes the norm $||\bm{\theta}||$ is likely the most
useful and is given by

$\bm{\theta}^{*}=\frac{{L}^{\dagger}\bm{\hat{\omega}}}{\sigma H^{\prime}(0)},$
(3)

where ${L}^{\dagger}$ is the Moore-Penrose pseudo-inverse of the Laplacian
matrixBen-Israel and Greville (1974). Importantly, to write down the exact
pseudo-inverse we require a full spectral decomposition of the Laplacian
matrix (i.e., global information of the network). For a possibly directed
network, the pseudo-inverse can be computed via the singular value
decomposition of $L$ as Skardal, Taylor, and Sun (2016)

${L}^{\dagger}=\sum_{n=2}^{N}\frac{\bm{v}_{n}\bm{u}_{n}^{T}}{\mu_{n}},$ (4)

where $0<\mu_{2}\leq...\leq\mu_{N}$ are the $N-1$ singular values of $L$ and
$\\{\bm{v}_{n}\\}$ and $\\{\bm{u}_{n}\\}$ are the set of right and left
singular vectors. Note that for the particular case of undirected networks the
singular values are given by the eigenvalues of $L$ and the left and right
singular vectors are are given by the eigenvectors of $L$, so $L^{\dagger}$ is
defined byt its eigenvalue decomposition.

Eq. (3) gives precisely the linearized synchronized state and therefore is
critical to understanding the interplay between structure and synchronization
dynamics in tasks involving the fully synchronized state, and for this reason
there is a myriad variety of phenomena and applications involving it. Eq. (3)
emerges in the derivation of objective functions used to optimize the degree
of synchronization in the systemSkardal, Taylor, and Sun (2014, 2016); Skardal
_et al._ (2016); Taylor, Skardal, and Sun (2016) and in problems related to
the control of desired synchronized states Skardal and Arenas (2015), critical
points and in network inference techniquesArenas, Díaz-Guilera, and Pérez-
Vicente (2006); Timme (2007); De Domenico (2017). Also, it is used to estimate
analytically the local and global susceptibility of a power-grid system
against structural failuresManik  _et al._ (2017) and to rank specific nodes
and edges according to its dynamical stabilityWitthaut and Timme (2012);
Coletta and Jacquod (2016). Furthermore, the object ${L}^{\dagger}$ provides,
by its own, a lot of insight on the network properties in diffusive processes
in generalVan Mieghem, Devriendt, and Cetinay (2017). In noisy dynamics, the
coupling of ${L}^{\dagger}$ with the covariance matrix of the input noise
determines the noise in the output signal such that the network can be
optimized to perform as a noise-canceling filterRonellenfitsch, Dunkel, and
Wilczek (2018). It plays a central role in the theory of electric circuits,
where it is used to obtain the resistance matrix that provides a notion of
diffusion distance between nodes and to determine the best spreader of the
networkVan Mieghem, Devriendt, and Cetinay (2017), as well as other node or
edge centrality measures, such as the random-walk betweennessNewman (2005).
From a mathematical point of view, ${L}^{\dagger}$ corresponds to the discrete
Green’s function of the networkChung and Yau (2000) since it quantifies how a
current injected in a node propagates to any other node of the networkEstrada
and Hatano (2008).

All the previous results are expressed in the literature in terms of the
spectra of $L$. The exact computation of ${L}^{\dagger}$ through the spectral
or singular value decomposition can be computationally costly and the
optimization and analysis of the interplay between structure and dynamics
relies on numerical schemes that are often treated as black-boxes. Therefore,
the networks studied in most of these cases are rather small and the
interesting properties that emerge in the optimal coupling of topology and
dynamics can only be observed _a posteriori_ , from the outcome of numerical
methods. Also, the computations and interpretations of the results require
global information of the network that is not available at the level of the
nodes, which usually operate in a decentralized manner Arenas  _et al._
(2008).

Motivated by the previous limitations, we derive a geometric expansion of Eq.
(3) in terms of increasingly further neighbourhoods of the nodes. This
approach allows for a decentralized, fast and accurate computation of
${L}^{\dagger}\bm{\omega}$, up to a desired degree of accuracy (or amount of
available information) and provides new analytical insights that allow for a
deeper understanding of the interplay between topology and dynamics in network
synchronization.

The remainder of this paper is organized as follows. In Sec. II we present our
main theoretical results, namely geometric series expansions for computing Eq.
(3), along with convergence proofs. In Sec. III we use local approximations,
generated from our theoretical results to identify link removals that improve
synchronization of networks. In Sec. IV we conclude with a discussion of our
results.

##  II Main Results

We begin by noting that rescaling the natural frequencies in Eq. (2) enables
us to set $\sigma$ and $H^{\prime}(0)$ both to one. To further simplify
notation we drop the hat-notation, in which case Eq. (2) reduces to Skardal,
Taylor, and Sun (2014); Skardal and Arenas (2015)

$\bm{\omega}=L\bm{\theta}.$ (5)

To solve Eq. (5) without the Moore-Penrose pseudo-inverse we focus our
attention more directly on $L$. First, using $L=D-A$, we write

$L=D(I-D^{-1}A).$ (6)

While $D$ is invertible (assuming that the network is connected and thus each
node has some positive degree), $(I-D^{1}A)$ is not. This can bee seen by
noting that $D^{-1}A$ is a stochastic matrix, and therefore has a leading
eigenvalue $\lambda_{1}=1$. Then $I-D^{-1}A$ has a zero eigenvalue, making it
singular. However, replacing $D^{-1}A$ with matrix $X$ that yields $(I-X)$
invertible, we have that

$[D(I-X)]^{-1}=(I-X)^{-1}D^{-1}.$ (7)

Moreover, the matrix $(I-X)^{-1}$ may be expanded in the geometric series

$(I-X)^{-1}=\sum_{m=0}^{\infty}X^{m}.$ (8)

The issue now arises that $X$ cannot be replaced by $D^{-1}A$, or more
specifically, we have that

$[D(I-D^{-1}A)]^{-1}\neq\sum_{m=0}^{\infty}(D^{-1}A)^{m}D^{-1},$ (9)

namely, on the left hand side the inverse is ill-posed, and this is reflected
by the fact that the series on the right-hand side diverges. However, this
does not rule out the possibility of the right-hand side converging when it is
applied to a vector of particular form. In fact, under the relatively mild
conditions of the network having a primitive adjacency matrix, when the series
is applied to an appropriately shifted frequency vector $\bm{\omega}$, the
right-hand-side does converge and yields a solution to Eq. (5), which leads to
the formulation of our first main result for undirected networks.

######  Theorem 1 (Convergence of the geometric series for undirected
networks).

Consider an undirected network with primitive adjacency matrix $A$ and a
frequency vector $\bm{\omega}$ with zero mean, i.e.,
$\langle\bm{\omega}\rangle=0$. Then, the infinite series

$\bm{\phi}=\sum_{m=0}^{\infty}({D}^{-1}{A})^{m}{D}^{-1}\bm{\omega},$ (10)

converges.

######  Proof.

We begin by denoting the symmetric normalized adjacency matrix as
$B=D^{-1/2}AD^{-1/2}$. Since $A$ is symmetric, so is $B$, and therefore its
normalized eigenvectors $\\{\bm{v}_{j}\\}_{j=1}^{N}$ form an orthonormal basis
for $\mathbb{R}^{N}$. Moreover, since $A$ is primitive, that is, there exists
some integer $M>0$ such that $A^{M}$ is strictly positive (note that this is
equivalent to $A$ being both irreducible and aperiodic), so is $B$. The Peron-
Frobenius theorem MacCluer (2000) then implies that $B$ has a single largest
eigenvalue $\lambda_{1}$ that is real and larger in magnitude than all other
eigenvalues, i.e., $\lambda_{1}>|\lambda_{j}|$ for $j=2,\dots,N$. Moreover,
since $B$ is normalized, we have $\lambda_{1}=1$ and $|\lambda_{j}|<1$ for
$j=2,\dots,N$. Finally, the leading eigenvector associated with the leading
eigenvalue $\lambda_{1}=1$ has entries that are proportional the the square
root of the degrees of the respective nodes, i.e.,
$\bm{v}_{1}\propto\bm{k}^{1/2}$.

Next, it is useful to define

$\bm{\phi}_{m}={D}^{-\frac{1}{2}}{B}^{m}{D}^{-\frac{1}{2}}\bm{\omega},$ (11)

so that the right-hand-side of Eq. (10) is given by
$\sum_{m=0}^{\infty}\bm{\phi}_{m}$. Defining the vector
$\bm{x}={D}^{-\frac{1}{2}}\bm{\omega}$, we now expand $\bm{x}$ via the
orthonormal basis of eigenvectors of $B$, namely,

$\bm{x}=\alpha_{1}\bm{k}^{\frac{1}{2}}+\alpha_{2}\bm{v}_{2}+...+\alpha_{N}\bm{v}_{N},$
(12)

where $\alpha_{i}=\langle\bm{x}_{i},\bm{v}_{i}\rangle$ are the coefficients
given by projections of $X$ onto the different eigenvector directions and we
assume that the eigenvector $\bm{v}_{1}=\bm{k}^{1/2}$ is also appropriately
normalized. Inserting Eq. (12) into Eq. (11) yields

$\begin{split}{D}^{\frac{1}{2}}\bm{\phi}_{m}&={B}^{m}\bm{x}\\\
&=\alpha_{1}\lambda_{1}^{m}\bm{k}^{\frac{1}{2}}+\alpha_{2}\lambda_{2}^{m}\bm{v}_{2}+...+\alpha_{n}\lambda_{n}^{m}\bm{v}_{n}.\end{split}$
(13)

Note now that for terms $j=2,\dots,N$, $\lambda_{j}^{m}$ decays geometrically
while $\lambda_{1}^{m}=1$. However, we now show that the coefficient
$\alpha_{1}$ must be zero. To see this, recall that the natural frequency
vector has mean zero, or in other words, $\bm{\omega}$ is orthogonal to the
constant vector $\bm{1}$, i.e., $\langle\bm{1},\bm{\omega}\rangle=0$. This is
equivalent to $\langle D^{1/2}\bm{1},D^{-1/2}\bm{\omega}\rangle$, or more
simply, $\langle\bm{k}^{1/2},\bm{x}\rangle=0$, which is precisely
$\alpha_{1}$. Thus, we have that

${D}^{\frac{1}{2}}\bm{\phi}_{m}=\alpha_{2}\lambda_{2}^{m}\bm{v}_{2}+...+\alpha_{n}\lambda_{n}^{m}\bm{v}_{n}.$
(14)

The convergence of the right-hand-side of Eq. (10) now follows quite easily:
since each of the finitely-many series
$\sum_{m=0}^{\infty}\alpha_{j}\lambda_{j}^{m}\bm{v}_{j}$ converges to
$\alpha_{j}\bm{v}_{j}/(1-\lambda_{j})$ for $j=2,\dots,N$, we have that the
full series converges to

$\begin{split}\phi&=D^{-1/2}\sum_{m=0}^{\infty}D^{1/2}\phi_{m}\\\
&=D^{-1/2}\sum_{m=0}^{\infty}\sum_{j=2}^{N}\alpha_{j}\lambda_{j}^{m}\bm{v}_{j}\\\
&=D^{-1/2}\sum_{j=2}^{N}\frac{\alpha_{j}\bm{v}_{j}}{1-\lambda_{j}},\end{split}$
(15)

which concludes the proof. ∎

A subtle issue that arises with the result given in Theorem 1 is that, while
the infinte series applied to the zero-mean vector $\bm{\omega}$ converges, it
does not necessarily converge to the minimum norm solution given by Eq. (3)
which turns out to the the zero-mean solution recovered by the Moore-Penrose
pseudo-inverse approach. This can be seen by noting that at initial truncation
of the geometric series we have that
$\langle\bm{\phi}_{0}\rangle=\langle\frac{\omega}{k}\rangle$, which is not
zero. However, this can be easily fixed by simply applying a constant shift to
the resulting vector to ensure that it is orthogonal to the constant
eigenvector $\bm{v}_{1}$ and the solution converges to the one with minimal
norm. In fact, this correction leads us to a generalization of our main
result, which importantly applies to both undirected and directed networks.

######  Theorem 2 (Convergence of the geometric series: General case).

Consider a network with primitive adjacency matrix $A$. Then, the infinite
series

$\displaystyle\bm{\phi}=\sum\limits_{m=0}^{\infty}\left(\bm{\phi}_{m}-\langle\bm{\phi}_{m}\rangle\right),$
(16)

where

$\displaystyle\bm{\phi}_{m}=(D^{-1}A)^{m}D^{-1}\bm{\omega},$ (17)

converges.

######  Proof.

The complete proof is shown in Appendix A.∎

Before proceeding to other practical questions and examples, we remark that
Theorem 2 includes no zero mean condition on the frequency vector nor is it
restricted to undirected networks, thus its potential usage is wider than
Theorem 1. In Fig. (1) we illustrate the utility of the more general case
given by Eqs. (16) and (17) when applied to a random network with random
allocated frequencies. We truncate the expansion at different neighborhood
orders $M$ and compare the approximated solution against the exact one given
by Eq.(3), showing that, for this particular random configuration, the
approximation is accurate even for small $M$.

Figure 1:  Scatter-plot of the exact Eq.(3) vs the truncated Eq.(17)
approximation of the stationary phases, for different neighborhood orders (M),
in a fixed Erdös-Rényi network of size $N=500$, mean degree $\langle
k\rangle=20$ with a normal distribution of frequencies $N(0,1)$ for three
different truncation orders.

Investigating our results described in Eq. (10) and Eqs. (16) and (17) more
deeply, we see that the geometric expansion expresses the solution of the
linear system as a sum of contributions of terms $\sim\omega/k$ coming from
increasingly further neighborhoods of the nodes. That is, the $m^{\text{th}}$
order term consists of terms $\omega_{j}/k_{j}$ corresponding to oscillators
located precisely $m$ links removed from a given oscillator. Thus, the series
is simply the Taylor expansion of the solution expanded about each oscillator
with higher-orders corresponding to larger network neighborhoods. It
constructs the exact solution by adding infinitely many incremental pieces of
local information in a polynomial basis that is not necessarily orthogonal.
Interestingly, this differs from the spectral decomposition used for
$L^{\dagger}$ in Eq. (3). In the latter, the solution is constructed by adding
$N-1$ pieces of global information (i.e. rank-one matrices for each non-zero
eigenvalue and its associated eigenvector in Eq.(4)), analogously to a Fourier
expansion that expresses the solution in the basis of orthogonal eigenmodes
with its associated eigenfrequencies (which carry global information on the
original function)McGraw and Menzinger (2007); Shuman  _et al._ (2013).

Furthermore, there is a clear connection between the applicability of this
method and Markov chains and random walks Masuda, Porter, and Lambiotte
(2017). In particular, for a Markov chain to have a unique, globally
attracting stationary state it is necessary and sufficient to have a primitive
transition matrix (or, equivalently, irreducible and aperiodic). For a network
with adjacency matrix $A$ the transition matrix for a random walk is given by
$D^{-1}A$ (or $AD^{-1}$, depending on the definition of the random walk). Here
$D^{-1}A$ appears in the geometric series because the solution is expressed in
terms of “in-neighborhoods”. Since this “spatial” expansion unfolds in terms
of in-neighborhoods of radius zero, one, two, three, etc., convergence of the
expansion implies that eventually these neighborhoods must include the whole
network, i.e., each node must be reachable from each node in a manner that is
eventually well-mixed, hence the need for the network adjacency matrix to be
irreducible and aperiodic, i.e., primitiveMasuda, Porter, and Lambiotte
(2017).

Turning now to a more computationally practical question, we investigate the
rate of convergence of this expansion. Assumign for simplicity the case of an
undirected network (so we may use Theorem 1), we consider the $M$-order
approximation defined as
$\bm{\theta}^{(M)}=\sum_{m=0}^{M}({D}^{-1}{A})^{m}{D}^{-1}\bm{\omega}$. Using
Eq. (13), the error is given by

$\begin{split}{D}^{\frac{1}{2}}(\bm{\theta}^{*}-\bm{\theta}^{(M)})&=\sum_{m=M+1}^{\infty}(\alpha_{2}\lambda_{2}^{m}\bm{v_{2}}+\cdots+\alpha_{n}\lambda_{n}^{m}\bm{v_{n}})\\\
&=\frac{\alpha_{2}\bm{v_{2}}}{1-\lambda_{2}}\lambda_{2}^{M+1}+\cdots+\frac{\alpha_{N}\bm{v_{N}}}{1-\lambda_{N}}\lambda_{N}^{M+1}.\end{split}$
(18)

As the order $M$ of the approximation increases, the dominant term in Eq. (18)
is that which corresponds to the second largest (in magnitude) eigenvalue
$\lambda_{2}$ of the normalized adjacency matrix $B$. (Recall that the largest
of eigenvalue of $B$ is $\lambda_{2}=1$.) Save for the unlikely scenario where
$\bm{x}=D^{-1/2}\bm{\omega}$ is exactly orthogonal to $\bm{v}_{2}$, in which
case $\alpha_{2}$ vanishes, then for large enough $M$ the mean square error
will scale geometrically with the magnitude of this second largest eigenvalue,
i.e.,

$||\bm{\theta}^{*}-\bm{\theta}^{(M)}||~\sim|\lambda_{2}|^{M+1}.$ (19)

Therefore, the smaller $|\lambda_{2}|$ is, the quicker the approximation will
converge. On the other hand, convergence will be slower for sparse networks
with either strong modularity or clustering as well as strong bipartite
structure, in which cases $\lambda_{2}$ tends to be close to $1$ and $-1$,
respectivelyNewman (2010); Masuda, Porter, and Lambiotte (2017) and a larger
number of terms are needed to obtain a desired level of accuracy. In Fig. (2)
we plot the error of the truncated approximation depending on the truncation
order for different values of the network density, modularity, clustering
coefficient and small-worldness and compare against the theoretical scaling
predicted by Eq.(19) observing a perfect agreement for sufficiently large
order $M$ in all the cases.

Figure 2:  Error of the approximation depending on the truncation order for
several interpolating network models, with fixed size $N=1000$. We study the
effect of density (varying the average degree), and fixing the mean degree at
$\langle k\rangle=10$, we also study the effect of modularity (interpolating
from a random network to a 2-modules configuration), clustering (using the
algorithm Holme and Kim (2002) to generate SF networks with controllable
clustering and the small-world model to interpolate between a ring and a
random network Watts and Strogatz (1998). In all the scenarios, we observe
that, for sufficiently high truncation order $M$, the error is dominated by
$\lambda_{2}$, the second largest eigenvalue of the normalized adjacency
matrix $B$. $\lambda_{2}$ increases with sparsity, modularity, and clustering
as expected.

##  III Local approximation of synchronization

Next we generate a local approximation from the geometric expansion presented
above to describe the degree of synchronization of a system given by the
Kuramoto order parameterKuramoto (2003b),
$re^{i\psi}=N^{-1}\sum_{j}e^{i\theta_{j}}$. First, after the linearizing about
a strongly synchronized state where all $|\theta_{i}|\ll 1$ (note that an
appropriate shift in initial conditions allows us to set the mean phase
$\psi=0$) one obtains Skardal, Taylor, and Sun (2014) $r\approx
1-\frac{1}{2\sigma^{2}N^{2}}||\bm{\theta}^{*}||^{2}$. Note that the degree of
synchronization increases, i.e., tends towards one, as the dispersion in the
phases is reduced. By truncating Eq. (10) at the first-order term and
neglecting the contribution of the shift to the mean, an approximation which
is accurate for locally tree-like networks for which $|\lambda_{j}|\ll 1$ for
$j=2,\dots,N$, we have that

$\theta_{i}\approx\frac{\omega_{i}}{k_{i}}+\left\langle\frac{\omega}{k}\right\rangle_{i}^{(1)}$
(20)

with
$\langle\frac{\omega}{k}\rangle_{i}^{(1)}=\frac{1}{k_{i}}\sum_{j=1}^{N}a_{ij}\frac{\omega_{j}}{k_{j}}$
being the average contribution of the first neighbours arriving at node $i$.
We can then directly write the order parameter as

$r\approx
1-\frac{1}{2\sigma^{2}N^{2}}\sum_{i=1}^{N}\left(\frac{\omega_{i}}{k_{i}}+\left\langle\frac{\omega}{k}\right\rangle_{i}^{(1)}\right)^{2}.$
(21)

The local unfolding of synchronization dynamics from the geometric expansion
that we use to write Eq. (21) allows us to gain analytical insight into the
interplay between topology and dynamics that improves synchronization as well
as understand several features that to date have only been investigated
numericallySkardal  _et al._ (2013); Skardal, Taylor, and Sun (2014, 2016);
Taylor, Skardal, and Sun (2016); Witthaut and Timme (2012); Coletta and
Jacquod (2016), for instance _i)_ linear degree-frequency correlations reduce
the absolute value of the phases at a given truncation order, which tend to
reduce the overall dispersion (increasing synchronization) Skardal  _et al._
(2013); Skardal, Taylor, and Sun (2016); Gómez-Gardeñes  _et al._ (2011),
_ii)_ negative frequency-frequency correlations between connected neighbors
tends to make the first order term of opposite sign (but smaller) to the local
term, and this reduces the dispersion (increasing synchronization) Skardal,
Taylor, and Sun (2014, 2016), _iii)_ synchronization is reduced when the
network become sparser. Even if the loss of links is compensated by an
increase of the weights, a sparse network will have more dispersion on the
phases with respect to a denser one, because in the latter the higher-order
terms sum over more neighbours and thus they vanish out more rapidly. This
provides a clear mechanistic interpretation of weight localization phenomena
Skardal, Taylor, and Sun (2019), and _iv)_ if frequencies are randomly
allocated, $r\sim 1-\langle k^{-2}\rangle$, where $\bm{k}$ are the strengths
(or degrees of the nodes). The inverse second moment is larger for
heterogeneous networks, thus homogeneity promotes synchronization in the
linearized regime Skardal, Taylor, and Sun (2019). This is opposite to what
occurs in the critical threshold, where $K_{c}\sim\langle k\rangle/\langle
k^{2}\rangle$, and homogeneity delays the critical threshold Arenas  _et al._
(2008). The interpretability of these effects emerge naturally from the local
description of synchrony, and although the particularities of each phenomenon
may require further analysis, the geometric unfolding allows to unveil the
underlying mechanistic rules that control the interplay of structure and
dynamics in synchronization.

Furthermore, the local approximation in Eq. (21) is useful in estimating the
impact that removing or adding a link has on the degree of synchronization. In
Ref. Taylor, Skardal, and Sun (2016) an approximation for the change in the
order parameter $r$ is given for undirected networks using a perturbation
analysis to approximate changes in the Laplacian eigenvalues, and in turn the
change to the order parameter. Using the geometric expansion this problem can
be addressed using only local information and directed networks may also be
treated. First, we denote the change in the order parameter $r$ as a result of
removing or adding a directed link coming from $q$ to $p$ as $\Delta r_{pq}$.
Assume that degrees are sufficiently large and no frequency or degree
correlations exist, it is straightforward to obtain from Eq. (21) that

$\Delta
r_{pq}\approx\mp\frac{1}{Nk_{p}}[\frac{\omega_{p}}{k_{p}}(\frac{\omega_{p}}{k_{p}}-\frac{\omega_{q}}{k_{q}})],$
(22)

where the $\mp$ sign accounts for removal/addition of links. First, the
absolute value $|\Delta r_{pq}|$ will be larger if node $p$, i.e., the node
with an incoming link from node $q$, has a large value of the ratio
$\omega_{p}/k_{p}$ (i.e., large frequency and/or low degree) and the
difference of phases (or ratios) between nodes $p$ and $q$ is large. Thus, the
most important links in terms of their contribution to the degree of
synchronization are those that connect nodes with large rations $|\omega/k|$.
Subsequently, connecting very similar nodes and hose that point towards nodes
with small $|\omega/k|$ contribute less to the degree of synchronization.

Figure 3:  Distribution of removed directed links $(p,q)$ in the phase-space
$(\omega_{p}/k_{p},\omega_{q}/k_{q})$ for an Erdös-Rényi network with $N=500$
and $\langle k\rangle=50$, and gaussian $g(\omega)$ with $N(0,1)$. The color
indicates whether if the removed link decreases (blue) or increases (red) the
degree of synchrony, calculated exactly from Eq.(3). Straight lines bound the
region predicted by Eq.(23) in the local approximation

Going one step further, using Eq. (22) we can predict which directed links
will produce the counter-intuitive effect of increasing (decreasing) synchrony
after its removal (addition). This effect, known as the Braess Paradox in the
context of road trafficNewman (2010), has been studied for oscillatory
networksWitthaut and Timme (2012); Coletta and Jacquod (2016); Motter and
Timme (2018) although the identification of these particular links relied on
numerical schemes or expressions in terms of the spectral decomposition of
$L$. Also, while these works study if the perturbations break the stability of
the current state, here we assume that the perturbations in the links drive
the system towards a new steady-state. In our formalism, we can directly
impose $\Delta r_{pq}>0$ for a link removal to obtain the condition

$\frac{\omega_{p}}{k_{p}}(\frac{\omega_{p}}{k_{p}}-\frac{\omega_{q}}{k_{q}})<0.$
(23)

The condition in Eq. (23) describes two regions of the
$(\omega_{p}/k_{p},\omega_{q}/k_{q})$-plane, namely the wedges
$\omega_{q}/k_{q}>\omega_{p}/k_{p}$ for $\omega_{p}>0$ and
$\omega_{q}/k_{q}<\omega_{p}/k_{p}$ for $\omega_{p}<0$. In terms of their area
these wedges describe a quarter of $(\omega_{p}/k_{p},\omega_{q}/k_{q})$
space. In Fig. (3) we plot the resulting change $\Delta r_{pq}$ for each
possible link removal in an Erdös-Rényi network of size $N=500$ with mean
degree $\langle k\rangle=50$, color-coding the change so that positive
(negative) changes are shaded more red (blue). We note that the positive
changes fit well within the wedges predicted by our local theory, which are
plotted in dashed black lines. In fact, approximately a quarter of directed
links have the potential to increase synchronization after its removal as
expected. While this phenomenon has been investigated in the context of
identical oscillatorsArenas  _et al._ (2008); Nishikawa and Motter (2006),
here the frequencies of the oscillators play a critical role in determining
which directed links are harmful or redundant, and we have shown that the
local approximation is sufficient to capture this phenomenon.

##  IV Conclusions

In this paper we introduced a geometric expansion that expresses the exact
solution of the linearized synchronization problem on a complex network as a
Taylor series in contributions from increasingly further network
neighborhoods. We proved that this series converges under relatively mild
conditions, namely, that the network has a primitive adjacency matrix, and
this method can be used for both directed and undirected networks. Moreover,
the approximation error from using a truncated series scales geometrically
with the second largest eigenvalue of the normalized adjacency matrix. We then
applied these theoretical results to gain analytical and mechanistic insight
on the optimal interplay of structure and dynamics on networks by deriving a
local approximation for the degree of synchronization in the system.
Furthermore, this approximation correctly identifies directed links may be
removed to improve synchronization properties, a phenomenon known as the
Braess paradox, as it captures the local interplay between network structure
and dynamics.

Our results are of both theoretical and practical relevance for the study of
synchronization on networks. The theory unfolds the stationary dynamics of the
system in the real (geometric) domain of the nodes and it provides a complete
understanding of the system from its local construction, unlike all previous
results found in literature that solve many specific problems in terms of
global spectral information. Also, we can use the theory to tackle important
practical problems such as network inferenceTimme (2007) and
optimizationSkardal, Taylor, and Sun (2014, 2016); Taylor, Skardal, and Sun
(2016), critical threshold estimationDörfler, Chertkov, and Bullo (2013),
control of synchronized states by local weight tuning or structure
modificationsSkardal and Arenas (2015); Arola-Fernández, Díaz-Guilera, and
Arenas (2018) and to predict how fluctuations in the weights or the
frequencies propagate to the macroscopic synchronization Arola-Fernández  _et
al._ (2020); Skardal, Taylor, and Sun (2019), just to name a few. Furthermore,
our theory can be applied beyond the synchronization problem, since it finds
utility in any forced diffusion system on networks that reduces to solving the
linear system of Eq. (5).

###### Acknowledgements.

LA-F and AA acknowledges the Spanish MINECO, Grant No. FIS2015-71582-C2-1.

## References

  * Prindle  _et al._ (2012) A. Prindle, P. Samayoa, I. Razinkov, T. Danino, L. S. Tsimring, and J. Hasty, “A sensing array of radically coupled genetic ‘biopixels’,” Nature 481, 39–44 (2012).
  * Witthaut and Timme (2012) D. Witthaut and M. Timme, “Braesss paradox in oscillator networks, desynchronization and power outage,” New J. Phys. 14, 083036 (2012).
  * Skardal, Taylor, and Sun (2014) P. S. Skardal, D. Taylor, and J. Sun, “Optimal Synchronization of Complex Networks,” Phys. Rev. Lett. 113, 144101 (2014).
  * Skardal, Taylor, and Sun (2016) P. S. Skardal, D. Taylor, and J. Sun, “Optimal synchronization of directed complex networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094807 (2016).
  * Taylor, Skardal, and Sun (2016) D. Taylor, P. S. Skardal, and J. Sun, “Synchronization of heterogeneous oscillators under network modifications: perturbation and optimization of the synchrony alignment function,” SIAM J. Appl. Math. 76, 1984–2008 (2016).
  * Skardal and Arenas (2015) P. S. Skardal and A. Arenas, “Control of coupled oscillator networks with application to microgrid technologies,” Sci. Adv. 1, e1500339 (2015).
  * Skardal and Arenas (2016) P. S. Skardal and A. Arenas, “On controlling networks of limit-cycle oscillators,” Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094812 (2016).
  * Dörfler, Chertkov, and Bullo (2013) F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks and smart grids,” Proceedings of the National Academy of Sciences 110, 2005–2010 (2013), https://www.pnas.org/content/110/6/2005.full.pdf .
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, “Collective dynamics of ’small-world’ networks,” Nature 393, 440–442 (1998).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Phys. Rep. 469, 93–153 (2008).
  * Pikovsky, Rosenblum, and Kurths (2003) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_ , 12 (Cambridge university press, 2003).
  * Kuramoto (2003a) Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * Newman (2010) M. Newman,  _Networks: An Introduction_ (Oxford University Press, Inc., New York, NY, USA, 2010).
  * Ben-Israel and Greville (1974) A. Ben-Israel and T. Greville,  _Generalized Inverses: Theory and Applications_ (Wiley, New York, 1974).
  * Skardal  _et al._ (2016) P. S. Skardal, D. Taylor, J. Sun, and A. Arenas, “Collective frequency variation in network synchronization and reverse pagerank,” Phys. Rev. E 93, 042314 (2016).
  * Arenas, Díaz-Guilera, and Pérez-Vicente (2006) A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, “Synchronization reveals topological scales in complex networks,” Phys. Rev. Lett. 96, 114102 (2006).
  * Timme (2007) M. Timme, “Revealing Network Connectivity from Response Dynamics,” Phys. Rev. Lett. 98, 224101 (2007).
  * De Domenico (2017) M. De Domenico, “Diffusion geometry unravels the emergence of functional clusters in collective phenomena,” Phys. Rev. Lett. 118, 168301 (2017).
  * Manik  _et al._ (2017) D. Manik, M. Rohden, X. Ronellenfitsch, H.and Zhang, S. Hallerberg, D. Witthaut, and M. Timme, “Network susceptibilities: Theory and applications,” Phys. Rev. E 95, 012319 (2017).
  * Coletta and Jacquod (2016) T. Coletta and P. Jacquod, “Linear stability and the Braess paradox in coupled-oscillator networks and electric power grids,” Phys. Rev. E 93, 032222 (2016).
  * Van Mieghem, Devriendt, and Cetinay (2017) P. Van Mieghem, K. Devriendt, and H. Cetinay, “Pseudoinverse of the Laplacian and best spreader node in a network,” Phys. Rev. E 96, 032311 (2017).
  * Ronellenfitsch, Dunkel, and Wilczek (2018) H. Ronellenfitsch, J. Dunkel, and M. Wilczek, “Optimal noise-canceling networks,” Phys. Rev. Lett. 121, 208301 (2018).
  * Newman (2005) M. E. J. Newman, “A measure of betweenness centrality based on random walks,” Social Networks 27, 39 – 54 (2005).
  * Chung and Yau (2000) F. Chung and S. T. Yau, “Discrete Green’s functions,” J. of Comb. Theory (A) 91, 191–214 (2000).
  * Estrada and Hatano (2008) E. Estrada and N. Hatano, “Communicability in complex networks,” Phys. Rev. E 77, 036111 (2008).
  * MacCluer (2000) C. R. MacCluer, “The many proofs and applications of perron’s theorem,” Siam Review 42, 487–498 (2000).
  * McGraw and Menzinger (2007) P. N. McGraw and M. Menzinger, “Analysis of nonlinear synchronization dynamics of oscillator networks by Laplacian spectral methods,” Phys. Rev. E 75, 027104 (2007).
  * Shuman  _et al._ (2013) D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE Signal Processing Magazine 30, 83–98 (2013).
  * Masuda, Porter, and Lambiotte (2017) N. Masuda, M. A. Porter, and R. Lambiotte, “Random walks and diffusion on networks,” Physics Reports 716-717, 1 – 58 (2017), random walks and diffusion on networks.
  * Holme and Kim (2002) P. Holme and B. Kim, “Growing scale-free networks with tunable clustering,” Physical review. E, Statistical, nonlinear, and soft matter physics 65, 026107 (2002).
  * Kuramoto (2003b) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Dover Publications, 2003).
  * Skardal  _et al._ (2013) P. S. Skardal, J. Sun, D. Taylor, and J. G. Restrepo, “Effects of degree-frequency correlations on network synchronization: Universality and full phase-locking,” Europhys. Lett.) 101, 20001 (2013).
  * Gómez-Gardeñes  _et al._ (2011) J. Gómez-Gardeñes, S. Gómez, A. Arenas, and Y. Moreno, “Explosive synchronization transitions in scale-free networks,” Phys. Rev. Lett. 106, 128701 (2011).
  * Skardal, Taylor, and Sun (2019) P. Skardal, D. Taylor, and J. Sun, “Synchronization of network-coupled oscillators with uncertain dynamics,” SIAM Journal on Applied Mathematics 79, 2409–2433 (2019).
  * Motter and Timme (2018) A. E. Motter and M. Timme, “Antagonistic Phenomena in Network Dynamics,” Ann. Rev. Cond. Matt. Phys. 9, 463–484 (2018).
  * Nishikawa and Motter (2006) T. Nishikawa and A. Motter, “Synchronization is optimal in nondiagonalizable networks,” Phys. Rev. E 73, 065106 (2006).
  * Arola-Fernández, Díaz-Guilera, and Arenas (2018) L. Arola-Fernández, A. Díaz-Guilera, and A. Arenas, “Synchronization invariance under network structural transformations,” Phys. Rev. E 97, 060301 (2018).
  * Arola-Fernández  _et al._ (2020) L. Arola-Fernández, G. Mosquera-Doñate, B. Steinegger, and A. Arenas, “Uncertainty propagation in complex networks: From noisy links to critical properties,” Chaos 30, 023129 (2020).

##  Appendix A Proof of Theorem 2

In this Appendix we now present the proof of Theorem 2.

######  Proof.

Rather than making use of the symmetric normalized adjacency matrix, as in the
Proof of Theorem 1, here we use the more classically stochastic matrix
$D^{-1}A$. However, this matrix shares similar properties with its symmetric
counterpart, namely, because $A$ is primitve, then so is $D^{-1}A$, and the
Peron-Frobenius theorem guarantees similar eigenvalue properties, namely there
is a single largest eigenvalue $\lambda_{1}$ that is real and larger in
magnitude than all other eigenvalues, i.e., $\lambda_{1}>|\lambda_{j}|$ for
$j=2,\dots,N$. Also, since $D^{-1}A$ is stochastic, $\lambda_{1}=1$ and
$|\lambda_{j}|<1$ for $j=2,\dots,N$. On the other hand, the leading
eigenvector is now given by the constant vector $\bm{v}_{1}\propto\bm{1}$.
Importantly, since $D^{-1}A$ is not symmetric (even if $A$ is), the
eigenvectors are not orthogonal to one another and cannot be used to form an
orthonormal basis for $\mathbb{R}^{N}$. Nonetheless, we may use these
eigenvalues as a (non-orthogonal) basis for $\mathbb{R}^{N}$ and uniquely
expand the vector $\bm{x}=D^{-1}\bm{\omega}$ using this basis, specifically

$\displaystyle\bm{x}=\sum_{j=1}^{N}\alpha_{j}\bm{v}_{j}.$ (24)

We then look at each eigenmode $j=1,\dots,N$ of the term $\bm{\phi}_{m}$,
namely $\bm{\phi}_{m}^{(j)}=\alpha_{j}(D^{-1}A)^{m}\bm{v}_{j}$ so that
$\bm{\phi}_{m}=\sum_{j=1}^{N}\bm{\phi}_{M}^{(j)}$. In terms of the full
expression for $\bm{\phi}$ [given in Eq. (16)] we then have that

$\begin{split}\bm{\phi}&=\sum\limits_{m=0}^{\infty}\sum_{j=1}^{N}\left(\bm{\phi}_{m}^{(j)}-\langle\bm{\phi}_{m}^{(j)}\rangle\right)\\\
&=\sum_{j=1}^{N}\sum\limits_{m=0}^{\infty}\left(\bm{\phi}_{m}^{(j)}-\langle\bm{\phi}_{m}^{(j)}\rangle\right).\end{split}$
(25)

We now treat the contribution of each eigenmode separately. We begin with the
eigenmodes $j\geq 2$ for which $|\lambda_{j}|<1$. First, we have that

$\displaystyle\bm{\phi}_{m}^{(j)}=\alpha_{j}(D^{-1}A)^{m}\bm{v}_{j}=\alpha_{j}\lambda_{j}^{m}\bm{v}_{j},$
(26)

and

$\displaystyle\langle\bm{\phi}_{m}^{(j)}\rangle=\langle\alpha_{j}\lambda_{j}^{m}\bm{v}_{j}\rangle=\alpha_{j}\lambda_{j}^{m}\langle\bm{v}_{j}\rangle,$
(27)

so together we have that

$\begin{split}\sum\limits_{m=0}^{\infty}\left(\bm{\phi}_{m}^{(j)}-\langle\bm{\phi}_{m}^{(j)}\rangle\right)&=\sum_{m=0}^{\infty}\alpha_{j}\lambda_{j}^{m}\left(\bm{v}_{j}-\langle\bm{v}_{j}\rangle\right)\\\
&=\frac{\alpha_{j}}{1-\lambda_{j}}\left(\bm{v}_{j}-\langle\bm{v}_{j}\rangle\right),\end{split}$
(28)

i.e., each component converges for $j\geq 2$.

To complete the proof we now show that the $j=1$ eigenmode, for which
$\lambda_{1}=1$, converges. In fact, it turns out that this component has no
contribution due to the shift of the mean. As in Eqs. (26) and (27), we have
that

$\displaystyle\bm{\phi}_{m}^{(1)}=\alpha_{1}(D^{-1}A)^{m}\bm{v}_{1}=\alpha_{1}\lambda_{1}^{m}\bm{v}_{1}=\alpha_{1}\bm{v}_{1},$
(29)

and

$\displaystyle\langle\bm{\phi}_{m}^{(1)}\rangle=\langle\alpha_{1}\lambda_{1}^{m}\bm{v}_{1}\rangle=\alpha_{1}\lambda_{1}^{m}\langle\bm{v}_{1}\rangle=\alpha_{1}\langle\bm{v}_{1}\rangle,$
(30)

so

$\displaystyle\bm{\phi}_{m}^{(1)}-\langle\phi_{m}^{(1)}\rangle=\alpha_{1}\left(\bm{v}_{1}-\langle\bm{v}_{1}\rangle\right),$
(31)

but since $\bm{v}_{1}\propto\bm{1}$, i.e., it’s constant, we have that
$\bm{v}_{1}=\langle\bm{v}_{1}\rangle$ and each term
$\phi_{1}^{(m)}-\langle\phi_{1}^{(m)}\rangle$ vanishes, which completes the
proof.∎