# The spark of synchronization in heterogeneous networks of chaotic maps

Antonio Montalbán  Department of Mathematics, University of California,
Berkeley, Berkeley 94720, USA  Rodrigo M. Corder<EMAIL_ADDRESS>Divisions of Epidemiology and Biostatistics, School of Public Health,
University of California, Berkeley, Berkeley 94720, USA

(August 28, 2024)

###### Abstract

We investigate the emergence of synchronization in heterogeneous networks of
chaotic maps. Our findings reveal that a small cluster of highly connected
maps is responsible for triggering the spark of synchronization. After the
spark, the synchronized cluster grows in size and progressively moves to less
connected maps, eventually reaching a cluster that may remain synchronized
over time. We explore how the shape of the network’s degree distribution
affects the onset of synchronization and derive an expression based on the
network’s construction that determines the expected time for a network to
synchronize. Understanding how the network structure affects the spark of
synchronization is particularly important for the control and design of more
robust systems that require some level of coherence between a subset of units
for better functioning. Numerical simulations in finite-sized networks are
consistent with this analysis.

###### pacs:

05.45.Xt, 89.75.Hc, 05.45.Ac

Synchronization stands as a pivotal phenomenon in networks, exerting a
profound impact across a spectrum of disciplines, including biology, chemistry
and physics, and various man-made systems. Notably, recent investigations
unveiled the phenomenon of chaotic units achieving sustained and stable
cluster synchronization within heterogeneous networks. Because many real-world
systems rely on cluster synchronization for their functioning, and natural
systems often comprise individuals with varying connection counts,
understanding the fundamental mechanisms that underlie the spark of synchrony
in heterogeneous networks holds promise for elucidating innovative control
strategies that amplify coherence among interacting entities. While extensive
research has addressed the conditions for synchronization, the dynamics
preceding the spark of synchrony remains elusive.

##  I Introduction

The complex structure of real-world networks has been extensively studied,
with significant progress made in recent years Milgram (1967); Barabási and
Albert (1999); Newman (2002); Radicchi et al. (2004); Zhou and Mondragón
(2004). These networks encompass a wide range of dynamical processes, from
natural systems in biology Barabasi and Oltvai (2004) to man-made
technological systems Pastor-Satorras et al. (2001). In particular, when a
network is composed of interacting dynamical systems, it can exhibit large-
scale coherent behavior that spontaneously emerges under certain conditions
Arenas et al. (2008). For instance, fireflies in a swarm are known to
synchronize their rhythms of flashing, resulting in highly correlated flashes
among a significant proportion of the swarm Sarfati et al. (2023). This
phenomenon, known as synchronization, has also been observed in various other
systems Winfree (2002).

Researchers have conducted extensive investigations into a class of models
called the Kuramoto models, which aim to describe the dynamics of coupled
oscillators Strogatz (2004); Rodrigues et al. (2016). Kuramoto himself
initially investigated the case of fully-connected networks with $N$ coupled
phase oscillators of equal strength Kuramoto (2013). His work revealed that,
in the continuum limit, there exists a critical coupling strength value -
dependent on the distribution of the phase oscillator frequencies - that
determines whether the phases of the oscillators in a coupled network will
remain incoherent or eventually evolve into synchronized behavior. Later, it
was found that the network topology has a significant influence on the
dynamics of these systems Lee (2005); Ichinomiya (2004); Gómez-Gardeñes et al.
(2007).

Extensive research has also been conducted on coupled collections of different
systems with more general dynamics, such as mixed chaotic and periodic
oscillators, chaotic maps, and others Baek and Ott (2004); Viana et al.
(2005); Pereira et al. (2017). In particular, in other work together with
Tiago Pereira and Zheng Bian Corder et al. (2023), we unravel the mechanism
for the emergence of cluster synchronization in heterogeneous random networks.
We developed a heterogeneous mean field approximation together with a self-
consistent theory to determine the onset and stability of the cluster. The
analysis showed that cluster synchronization occurs in a wide variety of
heterogeneous networks, node dynamics, and coupling functions. The system
dynamics before its asymptotic behavior, however, remains undisclosed in
random networks.

The transition to synchronization in coupled networks was studied from
different perspectives and it was observed that this phenomenon occurs or is
inhibited due to finite size fluctuations of the dynamical systems Komarov and
Pikovsky (2015); Ottino-Löffler and Strogatz (2018). In particular, it was
observed that the transition time to synchrony is exponentially distributed in
large homogeneous networks of chaotic circle maps Mendonca et al. (2023).
Here, we study the network dynamics before synchrony considering heterogeneous
networks of individuals whose dynamics are governed by Bernoulli maps Liang
and Cvitanović (2022). Notably, we do not assume any fixed form for the
network degree distribution and, for simplicity, we consider interactions
among coupled maps that can be described by a sinusoidal function Stankovski
et al. (2017). We depict fundamental processes that drive the spark of
spontaneous synchronization and observe that the transition time to synchrony
is also exponentially distributed in heterogeneous networks of chaotic maps.
We show that the theoretical analysis is consistent with numerical simulations
in finite-sized networks.

This work is organized as follows: in Section II, we define the dynamics of
the system, the network structure and some parameters that will be used to
measure the levels of synchronization in the network. Section III illustrates
the phenomenon under investigation with a computational simulation,
highlighting the crucial role of the connected oscillators in sparking
synchronization, and introduces the question we investigate through the paper.
We start discussing the phenomenon of the emergence of synchrony in Section IV
and, in Section V, we introduce the general model to describe the spark of
synchronization. In Section VI, we apply the model to describe the spark of
synchronization given a network structure. In Section VII, we present the
continuous approximation of the model, derive the formulas presented in
Section VI, and investigate how the shape of the degree distribution affects
the emergence of synchrony. In the last section we discuss our results.

##  II The model

We study networks with $N$ coupled maps $z_{i}$ satisfying,

$z^{t+1}_{i}\quad=\quad
2z^{t}_{i}+\frac{\alpha}{C}\sum_{j}A_{ij}\sin(z^{t}_{j}-z^{t}_{i})\quad\quad\mbox{mod
$2\pi$}.$ (1)

Here, $i=1,2,...,N$ labels each map, $\alpha$ is the network coupling
strength, $C$ is the network mean degree, $A_{ij}$ is 1 if nodes $i$ and $j$
are connected and 0 otherwise, and $z^{t}_{i}\in{\mathbb{R}}/2\pi{\mathbb{Z}}$
is the state of map $i$ at time $t$. We use $S^{1}$ to denote
${\mathbb{R}}/2\pi{\mathbb{Z}}$.

Let $\delta:{\mathbb{R}}^{+}\to{\mathbb{R}}^{+}$ be a probability density
function with mean 1. We consider random networks with degree distributions
$d\mapsto\delta(d/C)$. We denote by $d_{i}$ the degree of node $i$. In this
framework, the number of nodes $i$ with degree $d$ is $\delta(d/C)N/C$. Let
$w_{i}=d_{i}/C$. Thus, $w_{i}$ follow the distribution $\delta$ with mean 1.
Notice that, because we consider random graphs, the probability that node $i$
is connected to node $j$ is $d_{i}d_{j}/CN=w_{i}w_{j}C/N$.

We associate a complex number $u^{t}_{i}=e^{iz^{t}_{i}}$ to each state
$z^{t}_{i}$. In this paper, we will move interchangeably between a state
represented as a real number $z^{t}_{i}\in{\mathbb{R}}/2\pi{\mathbb{Z}}$ or a
complex number $u^{t}_{i}$ in the unit circle. For each map $i$, we define
$\vec{V}^{t}_{i}\in\mathbb{C}$, $r^{t}_{i}\in{\mathbb{R}}^{+}$ and
$\theta^{t}_{i}\in[0,2\pi]$ as follows,

$V^{t}_{i}=\sum_{j}A_{ij}u^{t}_{j}=r^{t}_{i}e^{i\theta^{t}_{i}}.$ (2)

By writing $\sin(z^{t}_{j}-z^{t}_{i})$ as the imaginary part of
$u^{t}_{j}\bar{u}^{i}_{j}$ (where $\bar{u}$ is the complex conjugate of $u$),
using Equation 2, we rewrite Equation 1 as,

$z^{t+1}_{i}=2z^{t}_{i}+\frac{\alpha}{C}\Im(V^{t}_{i}\bar{u}^{t}_{i}),$ (3)

where $\Im(v)$ denotes the imaginary part of the complex number $v$.

The distribution of the $V^{t}_{i}$’s. We approximate $V^{t}_{i}$ with a 2-D
Gaussian distribution for all maps $i$. Note that each $V^{t}_{i}$ is a sum of
$d_{i}$ vectors of the form $u^{t}_{j}$ (for $j$ connected to $i$), which are
distributed in a somewhat random manner around the circle. To calculate the
mean vector of those $u^{t}_{j}$’s, we sum all states of maps $j$, each
multiplied by the probability that map $j$ is connected to $i$, namely
$w_{i}w_{j}C/N$, and divide by the total number of maps $d_{i}=Cw_{i}$
connected to $i$. We obtain a vector that we call the weighted order parameter
and denote by $V^{t}$, i.e.,

$V^{t}=\frac{1}{N}\sum_{j}w_{j}\cdot u^{t}_{j}.$ (4)

Notice that $V^{t}$ is independent of $i$ and that the expected value for
$V^{t}_{i}$ is $d_{i}\cdot V^{t}$.

Next, to calculate the covariance matrix $\Sigma$ of $u^{t}_{j}$’s, we
consider the case when they are uniformly distributed in the circle (worst
case scenario). In this framework, on each coordinate, the variance of
$\sin(z)$ for $z$ uniformly in the interval $[0,2\pi]$ is $1/2$. The
covariance matrix $\Sigma$ is thus written as,

$\Sigma=1/2\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}.$

Therefore, by the central-limit theorem, we can estimate
$V^{t}_{i}\sim{\mathcal{N}}^{2}(d_{i}V^{t},(d_{i}/2){\bf I})$, where
${\mathcal{N}}^{2}$ is the Gaussian distribution on the plane and $\bf{I}$ is
the 2x2 identity matrix.

The weighted order parameter $V^{t}$. The weighted order parameter $V^{t}$ is
the average of $V^{t}_{i}$’s divided by the network mean degree $C$, i.e.,
$V^{t}=\frac{1}{CN}\sum_{i}\sum_{j}A_{ij}u^{t}_{i}$. It is computed by taking
into account the degree of each map $j$. Following the same procedure used to
define $V^{t}_{i}$, we can define $r^{t}\in{\mathbb{R}}^{+}$ and
$\theta^{t}\in[0,2\pi]$ as follows, $V^{t}=r^{t}e^{i\theta^{t}}$. Notice that
the most common approach in the literature is to compute the order parameter
as $\frac{1}{N}\sum_{j}u^{t}_{j}$, i.e., without taking into account the the
network degree distribution. In our heterogeneous framework, maps with higher
connectivity exert a greater influence on the network synchronization compared
to poorly connected maps. The phenomenon of synchronization in heterogeneous
networks is better captured by the weighted order parameter, as only a cluster
of maps with similar connectivity synchronizes at each time step. When such
clusters are small, little variation is observed in the unweighted order
parameter. In particular, this is more evident at the moment of spark, where
only a few highly connected maps spontaneously synchronize, making the
weighted order parameter more appropriate to capture this event.

##  III The phenomenon under investigation

Figure 1A presents a simulation exemplifying the phenomenon under
investigation in this work. We consider a heterogeneous network with
$N=100,000$ coupled maps with degree distribution $d$ randomly sampled from an
inverse-gamma distribution with mean $C=1,000$ and a power-law exponent
$\gamma=3$, and coupling strength $\alpha=15$. The initial states of the maps
are randomly sampled from a uniform distribution and the dynamics of the
system are governed by Equation 1. In this example, we observe the emergence
of spontaneous synchronization initially among few highly connected maps,
which occurs at a certain moment $t<10$. This moment marks the sparks of
synchronization. Notice in Figure 1A that the cluster of synchrony
progressively evolves over time, such that at each subsequent time step, the
most connected maps within the cluster lose synchronization concomitantly with
the synchronization of less connected maps. At a certain point, the system
reaches an asymptotic behavior where only maps with similar and intermediate
connectivity levels remain synchronized. Interestingly, throughout this
process, the maps within the cluster of synchrony exhibit a certain level of
coherence among themselves, but the cluster does not synchronize to a fixed
value. We refer to this phenomenon as partial synchronization—for a more
detailed discussion of the asymptotic behavior of this phenomenon, please see
Corder et al. (2023).

Figure 1: Emergency of partial synchronization in heterogeneous network of
coupled maps. To illustrate this phenomenon we consider a heterogeneous
network with $N=100,000$ coupled maps and degree distribution $d$ randomly
sampled from an inverse-gamma distribution with mean $C=1,000$ and a power-law
exponent $\gamma=3$ (truncated at 10,000), and coupling strength $\alpha=15$.
The figure presents a simulation where initial states of the maps
($z^{0}_{i}$, $i=1,2,...,N$) are randomly sampled from a uniform distribution
at the beginning of the simulation, and the dynamics of the system is governed
by Equation 1. Panel (A) presents the state of maps $z^{t}_{i}$ on the
$x$-axis and the normalized degree $w_{i}=d_{i}/C$ on the $y$-axis, with the
highest density areas in yellow. Panel (B) presents the level of
synchronization on the $x$-axis and normalized degree layer $w_{i}$ on the
$y$-axis, measured by the weighted order parameter $r^{t}_{w}$ (defined in
Section IV). We observe the emergence of spontaneous synchronization initially
among few highly connected maps at a certain moment $t<10$. The cluster of
synchrony progressively evolves over time, such that at each subsequent time
step, the most connected maps within the cluster lose synchronization
concomitantly with the synchronization of less connected maps. At a certain
point, the system reaches an asymptotic behavior where only maps with similar
and intermediate degree levels remain synchronized. Throughout this process,
maps within the cluster of synchrony exhibit certain level of coherence among
themselves, but the cluster does not synchronize to a fixed value. For a more
detailed discussion of the asymptotic behavior of this model, please see
Corder et al. (2023).

The question under investigation. Given a network design (defined by
parameters $N,C,\alpha,\delta$), what is the expected time for the spark to
occur when the network dynamics is governed by Equation 1? To answer this
question, we formulate a single number $s(N,C,\alpha,\delta)$ which
encapsulates all information about the network structure that allows to
estimate the expected time for the spark to occur. The number
$s(N,C,\alpha,\delta)$ is defined, as we will see in details in Section VI, as
follows:

$s(N,C,\alpha,\delta)=\frac{\alpha^{2}}{\sqrt{2N}}\ M^{\alpha}(\delta)\
K^{\alpha}_{\delta}(C),$

where
$M^{\alpha}(\delta)=M_{3}(\delta)\left(\sqrt{M_{2}(\delta)}+\frac{4}{\alpha\sqrt{M_{2}(\delta)}}\right)$,
with $M_{2}(\delta)$ and $M_{3}(\delta)$ representing, respectively, the
second and third cumulative moment of the probability density function
$\delta$. Note that $M^{\alpha}(\delta)$ depends solely on $\alpha$ and
$\delta$, and it measures the effect of the network degree distribution
sparsity on the expected time to spark. The term
$K^{\alpha}_{\delta}(C)=\frac{\int_{0}^{\infty}w^{3}e^{-w\frac{\alpha^{2}}{4C}}\delta(w)dw}{\int_{0}^{\infty}w^{3}\delta(w)dw}$
measures the impact of the noise resulting from the fact that each map is
connected to a distinct set of neighbors—the term $K^{\alpha}_{\delta}(C)$ is,
therefore, negligible for large values of $C$.

We will show in Section VI that the expected time to the spark of
synchronization is exponentially distributed and can be approximated by the
exponential of multiples of $s(N,C,\alpha,\delta)^{-2}$. An example of the
exponential behavior of the expected time to spark is presented in Figure 2A
where simulations were performed for different networks designs with degrees
randomly sampled from inverse-gamma distributions. The algorithm to calculate
the expected time to the emergence of partial synchrony in heterogeneous
networks of coupled maps is presented in the Appendix A.

Figure 2: Expected time for the emergence of partial synchronization in
heterogeneous network of coupled maps. Panel (A) presents the expected
sparking time of partial synchronization (in log-scale) as function of
$s(N,C,\alpha,\delta)^{-2}$. Each dot corresponds to the expected time for the
emergence of partial synchronization considering different network
configuration: $N\in\\{30000,50000,100000\\}$, $C\in\\{300,500,1000\\}$,
$\alpha\in\\{10,12,15\\}$ and $\gamma=\\{3,4,5\\}$ (power-law exponent of
inverse-gamma distribution $\delta$ with mean 1 and truncated at $10$—re-
scaled to keep the mean 1). The solid red line shows the expected expected
sparking time given by our theoretical results (${\mathbb{E}}(\sigma)$, as
defined in Section V). The log and inverse-square scales were chosen to
highlight the proportionality between coordinates, showing that the expected
time to spark behaves as an exponential of $s(N,C,\alpha,\delta)^{-2}$. The
dashed gray line shows the expected sparking time given by the theoretical
results with $\sigma$ re-scaled by 10%, (that is,
${\mathbb{E}}(1.1\cdot\sigma)$; more details in Section V). Panel (B) shows
the expected time to spark for a network with $N=30,000$ coupled maps and mean
degree $C=300$ computed for different values of $\alpha$ and $\gamma$. Dots
are obtained form network simulations and solid lines are the statistical
estimations (as described in Section V). Blue, orange, and green colors
represent $\gamma=3,4,5$, respectively.

##  IV The initiation of Synchronization

We stratify the process of synchronization in two parts:

  * •

The spark, that starts the process that is going to lead to synchronization,
and;

  * •

The build up, that occurs after the spark and leads to the asymptotic
synchronized behavior.

The spark. In Section II we defined the norm of the weighted order parameter
$V^{t}$, denoted as $r^{t}$, which quantifies the degree of order in the
system. The value of $r^{t}$ varies between 0 and 1, indicating complete
disorder or perfectly order, respectively, in the network dynamics. When the
system is out of synchrony, $r^{t}$ is close to zero and varies according to
an intrinsic noise inherent from the network construction, i.e., resulting
from the fact that few maps either synchronize or desynchronize purely by
chance before the emergence of synchronization. After the spark of
synchronization, $r^{t}$ might progressively increase until it reaches an
asymptotic behavior. It is important to note that the sparking of
synchronization is not a well-defined moment. However, in our simulations, we
observe that there are sufficiently large thresholds, such that when $r^{t}$
exceeds that threshold, the system is highly likely to evolve towards
synchronization. The ”spark” happens at some point before $r^{t}$ reaches the
threshold, and it can be defined and studied in statistical terms, as we will
explain in more detail below. The ”spark” is the phenomenum discussed
throughout all sections of this work.

The build up. A system may be out of synchrony configuration for many time
steps before a moment in which, spontaneously, maps start transitioning to a
more synchronized dynamics. After this moment, the level of synchronization
starts increasing for a few time steps before reaching the asymptotic
behavior. To understand this process, first, we stratify maps by their
connectivity level $w$ and compute the associated $r^{t}_{w}$ (norm of the
weighted order parameter of maps with connectivity $w$). We notice that at the
beginning of the transition phase from chaos to partial synchrony, the larger
contribution to $r^{t}$ comes from nodes of higher degree (Figure 1B). As the
system evolves, the cluster of synchrony progressively moves through less
connected oscillators until reaching maps that may remain synchronized over
time. This phenomenon is due to the fact that, at each time step $t$, maps
with connectivity $w$ close to the ratio $2/\alpha r^{t}$ are more likely to
synchronize (more details about this condition is presented in the Appendix
B). If they synchronize, the level of synchronization in the system increases,
i.e., $r^{t}$ increases. Thus, in the next time step, maps with lower
connectivity $w$ will be more likely to synchronize since the fraction
$2/\alpha r^{t}$ is reduced. As $r^{t}$ keeps increasing, the level $w$ at
which maps are more likely to synchronize keeps decreasing. At a certain
moment, the cluster of synchronization comprises maps with a certain
connectivity level $w$ in which $r^{t}$ is close to its maximum value for that
system. At that point, $r^{t}$ stops increasing and maybe even decrease, and
the system is close to its asymptotic behavior. We analyzed the asymptotic
properties of this phenomenon in details together with Tiago Pereira and Zheng
Bian in another work Corder et al. (2023).

##  V The Markov process

The main claim of this paper is that the dynamics of the weighted order
parameter $V^{t}$ prior to partial synchronization can be closely approximated
by a Markov process as in the Definition 1.

######  Definition 1.

Given parameters $k,\sigma\in{\mathbb{R}}^{+}$, let

  * •

$v_{0}=0\in{\mathbb{C}}$ and

  * •

$v_{t+1}=kv_{t}^{2}+\epsilon_{s}$,

where $\epsilon_{s}\in{\mathbb{C}}$ follows a 2D-Gaussian distribution with
mean 0 and covariance matrix $\sigma^{2}{\bf I}$, where ${\bf I}$ is the 2x2
identity matrix. Here, the square operation on $v_{t}$ refers to the square
operation on complex numbers.

Notice that, if $||v_{t}||\ll 1/k$, then $||v_{t+1}||\approx\epsilon_{s}$ for
large enough Gaussian-noise $\epsilon_{s}$. In this scenario, the dynamics of
$v_{t}$ will oscillate around the origin according to the Gaussian-noise.
Eventually, $||v_{t}||$ may increase and the term $kv_{t}^{2}$ will have more
influence on next $v_{t+1}$. If at some point $||v_{t}||$ is much greater than
both $1/k$ and $\sigma$, then $k||v_{t}^{2}||$ will be even greater. After
this moment, $v_{t}$ will very likely keep increasing forever and escape any
bounded region. We are interested in studying how long it takes for this
escape to happen.

An important observation is that the Markov process under investigation is
determined only by the product $k\sigma$. Suppose a second process governed by
parameters $k^{\prime}$ and $\sigma^{\prime}$ such that
$k^{\prime}\sigma^{\prime}=k\sigma$. It is not hard to see that the dynamics
of the two systems are equivalent via the transformation
$v^{\prime}=(k/k^{\prime})v$, as the noise in the second system would have
standard deviation $(k/k^{\prime})\sigma=\sigma^{\prime}$. For the rest of
this section we consider the case for $k=1$.

Figure 3 shows that the expected time to escape for the Markov process as in
Definition 1, with a fixed value of $\sigma$ (and $k=1$), presents a similar
pattern for the expected time to escape obtained from the dynamics of the
network system. We proceeded as follows. Let’s consider a value $B$ large
enough such that, if the $||v_{t}||>B$, we can be almost sure that the process
will never return to a neighborhood of 0 (say $B=100$). Consider that the
random variable $T$ models the time $t$ in which $||v_{t}||>B$. Notice that,
by time $T$, the process has already escaped the neighborhoods of 0 (almost
surely) forever. Notice also that the event of escaping happens before $T$,
though we do not exactly know the precise moment. However, besides some noise,
we obtain from simulations that the distribution of $T$ is very close to a
constant $c$ plus an geometric distribution. Importantly, the decay $\lambda$
of the geometric distribution is independent of $B$. The value of the constant
$c$ is given by the number of time steps from the escaping moment itself to
the moment in which $||v_{t}||>B$. Therefore, the constant $c$ depends on the
choice of $B$. Figure 3 presents an example with $B=100$ (and $k=1$) which
graphically implies $c=5.25$ and $\lambda=12.01$.

Figure 3: Expected time of escaping (${\mathbb{E}}(\sigma)$). Panels (A) and
(C) show the the probability mass function and the expected value of $T-x$
given that $T>x$ plotted against $x$, where $T$ is the time to escape random-
variable for the Markov process (Def. 1) with $\sigma=0.36$, $B=100$. The
orange and green dashed lines represent $E(T_{g}-x|T_{g}>x)$ and
$E(T_{c}-x|T_{c}>x)$ where $T_{c}=17.26$ and $T_{g}$ is a geometric random
variable with mean $\lambda=12.01$. This shows how $T$ is very close to
$5.25+T_{g}$. Panel (B) and (D) show the probability mass function and the
same expected value of $T-x$, but with $T$ representing the time-to-escape
obtained from a simulation of a network with $N=50,000$, $C=500$, $\gamma=4$
and $\alpha=12$, and dynamics governed by Equation 1. Orange and green dashed
lines consider values $T_{c}=15.41$ and $\lambda=11.42$.

To show that the distribution of the random variable $T$ is close to a
constant plus a geometric distribution, in Figure 3 we illustrate the expected
value $E(T-x|T>x)$ of the remaining number of time steps, given that the
number of steps is above $x$. One can see that from some fixed value of $x$
onward, $E(T-x|T>x)$ is constant and, thus, $T$ is memoryless like the
geometric distribution. Recall that the geometric is the only distribution on
$\mathbb{N}$ for which $E(T-x|T>x)$ is independent of $x$. Notice that, for
small values of $x$, $E(T-x|T>x)$ is linear in $x$ with slope $-1$.

######  Definition 2.

We define the expected time to escape as the expected value of the geometric
part of the distribution of a random variable $T$, i.e., the value of
$E(T-x|T>x)$ for some large enough $x$. We call this function
${\mathbb{E}}(\sigma)$.

Importantly, we do not derive a closed form equation for the function
${\mathbb{E}}(\sigma)$.

##  VI The main result

Here is our main claim, but now accounting for the network structure and its
dynamics: Given a network design defined by parameters $N$, $C$, $\alpha$ and
$\delta$, the dynamics of the weighted order parameter $V^{t}$ before
synchronization behaves as the Markov process as in the Definition 1 with
parameters

$\begin{split}\sigma&=\sqrt{M_{2}(\delta)/(2N)}\text{, and}\\\ k&=\alpha^{2}\
K^{\alpha}_{\delta}(C)\ M_{3}(\delta)(1+4/\alpha M_{2}),\end{split}$ (5)

where $M_{2}(\delta)=\int_{0}^{\infty}w^{2}d\delta(w)$ and
$M_{3}(\delta)=\int_{0}^{\infty}w^{3}d\delta(w)$ are as defined before, the
second and third cumulative moment of the probability density function
$\delta$, respectively. $K^{\alpha}_{\delta}(C)$ (derived in Section VII) is
the only parameter that depends on $C$ and it measures the impact of the noise
resulting from the fact that each map is connected to a distinct set of
neighbors. As an example, Figure 2A presents simulations performed for
different networks designs with degrees randomly sampled from inverse-gamma
distributions and shows that the expected time to synchronization can be
approximated by the Markov process as in the Definition 1, with parameters
defined based on the network structure as presented in 5.

The expected time to spark of synchronization can thus be approximated by
${\mathbb{E}}(k\sigma)$, which is approximated by

${\mathbb{E}}\left(\frac{\alpha^{2}}{\sqrt{2N}}\ M^{\alpha}(\delta)\
K^{\alpha}_{\delta}(C)\right),$ (6)

where $M^{\alpha}(\delta)$ (derived in Section VII), is a term that measures
the effect of the sparsity of the network degree distribution $\delta$ on the
expected time to spark. Figure 4 presents how $M^{\alpha}(\delta)$ and
$K^{\alpha}_{\delta}(C)$ vary according to $\delta$ and $C$.

Figure 4: The parameters of the Markov process under investigation. Panels (A)
and (B) present the shape of inverse gamma distributions in linear-linear and
log-log scales, respectively, for power-law exponents $\gamma=3$ (blue curve),
$\gamma=4$ (yellow) and $\gamma=5$ (green). Panel (C) presents how
$M^{\alpha}(\delta)$ varies according to the power-law exponent $\gamma$.
Panel (D) presents how $K^{\alpha}_{\delta}(C)$ varies according to $C$. In
panels (C) and (D), blue, orange and green curves correspond to $\alpha=10,12$
and $15$.

Equation 6 shows how the values $\alpha$, $C$, $N$, and the probability
distribution $\delta$ are related to the expected time to spark. Some
observations: (i) the larger the value of the coupling strength $\alpha$, the
faster the spark; (ii) the larger the number of maps $N$ in a network, the
slower the spark; (iii) the time to spark is correlated with the square root
of the number of maps $N$ and to the square of the network coupling strength
$\alpha$; (iv) the more disperse the network degree distribution (larger
values of $M_{3}(\delta)$), the faster the spark; (v) the term
$\sqrt{M_{2}(\delta)}+\frac{4}{\alpha\sqrt{M_{2}(\delta)}}$ is smallest when
$M_{2}(\delta)=4/\alpha$. Since $M_{2}(\delta)>1>4/\alpha$, this term also
increases with $M_{2}(\delta)$ (dispersion of $\delta$).

Notice that the parameter $\sigma$ is the standard deviation of the weighted
order parameter $V^{t}$, calculated for $N$ maps with states randomly sampled
from a uniform distribution, and on a network with degree distribution
$\delta$. Therefore, $\sigma=\sqrt{M_{2}(\delta)/(2N)}$. We argue that this is
a good approximation since the state of maps before synchrony can be
approximated by a uniform distribution. In the next section we deeply
investigate the dynamics of the system and present the derivation of the
formula for $k$.

##  VII The continuous approximation

In this section, we dive into a deep investigation of the network dynamics of
coupled maps governed by Equation 3 and elucidate the motivations behind the
claims made in previous sections. We will derive a continuous approximation to
Equation 3 such that the states of maps are well approximated by a probability
distribution that evolves over time. We notice, however, that we are not
deriving the limit as the size $N$ of the graphs converges to $\infty$.
Instead, we assume that $N$ is sufficiently large such that the continuous
version of the model provides a good approximation of its dynamics, but not so
large as to eliminate the stochastic behavior inherent in a finite-sized
network. We will demonstrate that small random variances in the degree
distribution introduced by random networks play a crucial role in the sparking
of synchronization.

Acting on continuous probability distributions. For each connectivity layer
$w\in{\mathbb{R}}^{+}$, we consider a probability density function
$\rho_{w}(z)$ with $z\in S^{1}$. We use $\rho$ to denote the probability
density for the pair $(w,z)\in{\mathbb{R}}^{+}\times S^{1}$, denoted as
$\rho(w,z)=\rho_{w}(z)\delta(w)$. This function evolves on time (we omit the
superscript $t$) as the dynamical system acts on the space of probability
density functions $\rho$ on ${\mathbb{R}}^{+}\times S^{1}$ that are compatible
with the network degree distribution $\delta$.

Consider a fixed value of $w$ and denote $\tau(z)=\rho_{w}(z)$ as a
probability density function on $S^{1}$. Given the order parameter
$V\in{\mathbb{C}}$, we define an operator $\mathbb{F}_{V}$ acting on
probability distribution functions defined on $S^{1}$. $\mathbb{F}_{V}(\tau)$
describes the probability distribution on $S^{1}$ obtained by applying one
time step of the dynamics. To this end, for each map $i$ with connectivity $w$
and state $z_{i}\in S^{1}$, we apply the function,

$z_{i}\mapsto 2z_{i}+\frac{\alpha}{C}\Im(V_{i}\bar{u}_{i}),$

where $\bar{u}_{i}=e^{-iz_{i}}$ and $V_{i}$ is the sum of all states
$u_{j}$’s, with $j$ connected to $i$ (as described in Section II). Recall from
Section II that $V_{i}$ is distributed in the complex plane according to a 2D
Gaussian distribution with mean $d_{i}V$ and covariance matrix $(d_{i}/2){\bf
I}$, where $d_{i}=wC$ is the number of connections of map $i$. It implies that
$\Im(V_{i}\bar{u}_{i})$ is distributed following a 1D Gaussian distribution
with mean $wC\Im(V\bar{u}_{i})$ and variance $wC/2$. Thus, the dynamics can be
described by mapping $z_{i}\mapsto 2z_{i}+\alpha w\
\Im(V\bar{u}_{i})+\epsilon_{i}$, where $\epsilon_{i}$ is a 1D Gaussian random
variable with mean $0$ and variance $\alpha^{2}w/2C$.

We decompose the operator $\mathbb{F}_{V}$ in two operators also acting on
probability distribution functions defined on $S^{1}$:

$\mathbb{F}_{V}=\mathbb{D}_{\sqrt{\alpha^{2}w/2C}}\circ{\mathbb{M}}_{\alpha
wV}.$

The main operator ${\mathbb{M}}_{P}$, for $P\in{\mathbb{C}}$, maps a
distribution $\tau$ to the distribution obtained after applying the function
$z_{i}\mapsto 2z_{i}+\Im(P\bar{u}_{i})$. The diffusion operator
$\mathbb{D}_{\nu}$ maps a distribution $\tau$ to the distribution we get from
applying the function $z_{i}\mapsto z_{i}+\epsilon_{i}$, where $\epsilon_{i}$
is a normal random variable with mean $0$ and standard deviation $\nu$. In
this framework, the mean field $E_{\tau}$ at each layer $w$ can be obtained,
at each time step $t$, as in Def. 3.

######  Definition 3.

Given a probability distribution function $\tau$ on $S^{1}$, let $E_{\tau}$ be
the expected value of $e^{iz}$ where $z\sim\tau$, that is,

$E_{\tau}=\int_{S^{1}}e^{iz}\tau(z)dz.$

Next, we analyze the operators $\mathbb{D}_{\nu}$ and ${\mathbb{M}}_{P}$.

The diffusion operator $\mathbb{D}_{\nu}$. $\mathbb{D}_{\nu}(\tau)$ is the
distribution obtained after mapping a distribution $\tau$ to the function
$z_{i}\mapsto z_{i}+\epsilon_{i}$, where $\epsilon_{i}$ is a 1D Gaussian
random variable with mean $0$ and standard deviation $\nu$. The diffusion
operator acts on the mean field $E_{\tau}$ by shrinking its length by a factor
of $e^{-\nu^{2}/2}$, but keeping its direction (Lemma 4).

######  Lemma 4.

$E_{\mathbb{D}_{\nu}(\tau)}=e^{-\nu^{2}/2}E_{\tau}$ (proof in the Appendix C).

The main operator ${\mathbb{M}}_{P}$. We study the main operator
${\mathbb{M}}_{P}$ on probability distributions that are very close to the
uniform distribution, i.e., the case before the spark of synchronization.
Since $z\in[0,2\pi]$, the uniform distribution is given by
$\tau_{0}(z)=1/2\pi$.

######  Definition 5.

Given $a=\epsilon e^{i\theta}\in{\mathbb{C}}$, let $\tau_{a}$ be the
probability distribution on $S^{1}$ of the form

$\tau_{a}(z)=\frac{1}{2\pi}\big{(}1+2\epsilon\cos(z-\theta)\big{)}.$

Notice that we can also write
$\tau_{a}(z)=\frac{1}{2\pi}\big{(}1+2\Re(e^{iz}\bar{a})\big{)}$, where
$\bar{a}$ is the complex conjugate of $a$ and thus the uniform distribution
corresponds to $\tau_{a}$ for $a=0$. Interestingly, we will see below that the
application of ${\mathbb{M}}_{P}$ or $\mathbb{D}_{\nu}$ on a probability
distribution of the form $\tau_{a}(z)$ maps it to another distribution of the
form $\tau_{a^{\prime}}(z)$ up to a small error. Before, we show on Lemma 6
that the expected value of $\tau_{a}$ is $a$.

######  Lemma 6.

$E_{\tau_{a}}=a$ (proof in the Appendix C).

The distribution $\tau_{a}$ is essentially the simplest distribution on
$S^{1}$ that one can think of with mean field $a$. Thus, we assume that,
before synchronization, there exist a complex number $a_{w}$ for each $w$,
with norm of $a_{w}$ small, such that $\rho_{w}=\tau_{a_{w}}$. Lemmas 7 and 8
present how the operators ${\mathbb{M}}_{P}$ and $\mathbb{D}_{\nu}$,
respectively, act on distributions of the form $\tau_{a}$.

######  Lemma 7.

Given $a\in{\mathbb{C}}$ and $P\in{\mathbb{C}}$ ,

${\mathbb{M}}_{P}(\tau_{a})\approx\tau_{a^{\prime}},$

where $a^{\prime}=P(P+4a)/8$ (proof in the Appendix C).

######  Lemma 8.

Given $a\in{\mathbb{C}}$ and $\nu\in{\mathbb{R}}^{+}$,

$\mathbb{D}_{\nu}(\tau_{a})\approx\tau_{a^{\prime}},$

where $a^{\prime}=e^{-\nu^{2}/2}a$ (proof in the Appendix C).

The dynamics. Here we consider all connectivity layers $w$ and compute the
dynamics of the order parameter $V$. We can calculate the order parameter at
the next time step, $V^{+}$, and show that the map $V\mapsto V^{+}$ can be
well approximated by the Markov process defined in Section V. Since we are
dealing with the dynamics of the network before the sparking of
synchronization, we assume that all distributions $\rho_{w}$ are close to
uniform. For each $w$, let $a_{w}=E_{\rho_{w}}$ (close to $0$ before
synchronization). Notice that $a_{w}$ must satisfy,

$\int wa_{w}\delta(w)dw=V.$

We also assume that $\rho_{w}$ is of the form $\tau_{a_{w}}$, for all $w$, and
calculate $\rho^{+}_{w}$ (density $\rho_{w}$ at the next time step) as
follows,

$\rho^{+}_{w}=\mathbb{F}_{V}(\rho_{w})=\mathbb{D}_{\alpha\sqrt{w/2C}}({\mathbb{M}}_{\alpha
wV}(\tau_{a_{w}}))=\tau_{a^{\prime}_{w}}$

where $a^{\prime}_{w}=e^{-\alpha^{2}w/4C}\alpha wV(\alpha wV+4a_{w})/8$.
Therefore, the expected value for $V^{+}$ is obtained by,

$\displaystyle V^{+}$ $\displaystyle=$ $\displaystyle\int
w\left(e^{-\alpha^{2}w/4C}\alpha wV(\alpha wV+4a_{w})/8\right)\delta(w)dw$ (7)
$\displaystyle=$
$\displaystyle\alpha^{2}(V)^{2}\int\left(w^{3}e^{-\alpha^{2}w/4C}\delta(w)/8\right)dw$
$\displaystyle+\alpha
V\int\left(w^{2}a_{w}e^{-\alpha^{2}w/4C}\delta(w)/2\right)dw.$

Based on Equation 7, we define $K^{\alpha}_{\delta}(C)$ as,

$K^{\alpha}_{\delta}(C)=\frac{\int_{0}^{\infty}w^{3}e^{-w\frac{\alpha^{2}}{4C}}\delta(w)dw}{\int_{0}^{\infty}w^{3}\delta(w)dw}.$

This term is relevant only for small values of $C$. If the mean degree $C$ is
large, $e^{-w\frac{\alpha^{2}}{4C}}\approx 1$ and thus
$K^{\alpha}_{\delta}(C)\approx 1$. For small values of $C$,
$K^{\alpha}_{\delta}(C)$ is always less than $1$, i.e., it slows down the
sparking time. The reason for this lies in the fact that each map $i$ receives
signals from a distinct set of $w_{i}C$ neighbors, resulting in some random
noise in the received signals for each map. This noise has a detrimental
effect on the synchronization of the system.

The second term of Equation 7 presents the random variable $a_{w}$. We show in
the Appendix D that $\int w^{2}a_{w}e^{-\alpha^{2}w/4C}\delta(w)$ has expected
value $VK^{\alpha}_{\delta}(C)/M_{2}$.

The variance on the order parameter. After computing the distribution
$\rho^{+}$, we take $N$ random points $(w_{i},z_{i})$. The next order
parameter is then obtained as $V^{+}=\sum_{i=1}^{N}w_{i}e^{iz_{i}}$. The
expected value of $w_{i}e^{iz_{i}}$ is
$V^{\prime}=(\alpha^{2}K^{\alpha}_{\delta}(C)V^{2}/8)(1+4/\alpha M_{2})$. As
for the variance, we obtain $E(w_{i}^{2})Var(e^{iz_{i}})$. Since $\rho^{+}$ is
very close to the uniform distribution, we can approximate
$Var_{\rho^{+}}(e^{iz_{i}})\approx Var_{unif}(e^{iz_{i}})=(1/2){\bf I}$, where
${\bf I}$ is the 2x2 identity matrix. Moreover,

$E(w_{i}^{2})=\int w^{2}\delta(w)dw=M_{2}(\delta)$

is the second momentum of the degree distribution $\delta$. Taking the mean of
$N$ points, we can use the central-limit theorem and conclude that $V^{+}$ is
2D-Gaussian distributed in the complex plane with mean $V^{\prime}$ and
variance $M_{2}(\delta)/2N$.

##  VIII Discussion

We studied the mechanisms for the spark of synchronization in heterogeneous
networks of Bernoulli coupled maps. We obtained a continuum limit
approximation of the model and derived an expression whose expected value
determines the sparking time of synchrony given a network construction. We
showed with numerical simulations in finite-sized networks that our results
are consistent with the analysis.

The correlation of the expected sparking time of synchrony and the network
construction (Equation 6) shows that, for a given coupling strength, mean
degree and degree distribution, the larger the network, the longer is the
expected time the system takes to synchronize. Furthermore, the time it takes
to spark is more influenced by changes in the coupling strength (powered-two
correlation) rather than the number of maps in the network (squared-root
correlation). One can also see that the more heterogeneous the network degree
distribution (i.e., the larger $M_{2}(\delta)$ and $M_{3}(\delta)$), the
sooner the system is expected to synchronize. Also, larger mean degrees lead
to faster synchronization since it implies larger values of
$K^{\alpha}_{\delta}(C)$; however, the increase in $K^{\alpha}_{\delta}(C)$ is
rather slow for larger values of $C$ and its variation after a certain
threshold is expected to have small effects on the time to synchronization.

Our results are consistent with others Pereira et al. (2017, 2017); Restrepo
et al. (2006) and shows that the emergence of coherence in large collections
of heterogeneous coupled chaotic systems is also highly dependent on the
network topology. However, while most analysis of the synchronization in
coupled networks with interacting dynamical systems investigates the critical
coupling strength, we looked at the problem from a new perspective. By
studying the dynamics of systems on heterogeneous networks before synchrony,
we derived an expression that correlates the parameters of the network
structure with the expected time to synchronization. It is important to note
that due to the intrinsic error inherent in finite-sized networks being
modeled with normal distributions, there is always a chance that a network may
lose synchronization. However, the network dynamics may spontaneously return
to synchrony quickly for some identifiable configurations. This analysis could
be pertinent for various applications where synchrony among a set of
individuals is a desirable condition, particularly those where one can control
the network’s structure.

To our knowledge, this is the first attempt to investigate the expected time
for a network to reach synchrony considering heterogeneous networks without a
fixed form for the degree distribution. Our simulations show that the Markov
process presented in this work is a suitable tool to model the sparking of
synchronization in heterogeneous networks of coupled maps. Our results,
however, present a small deviation (approximately 10%; Figure 2A) which does
not affect qualitatively the conclusions derived from the analysis. Future
investigations may try to understand the reasons for this deviation as well as
derive a closed-form equation for the expected time to synchronization.
Whether to study the dynamics of neurons Iaccarino et al. (2016); Martorell et
al. (2019), power grids Motter et al. (2013) or any other heterogeneous
networks of coupled units, we argue that understanding the fundamental
mechanisms for the spark of synchrony can shed light into new control
strategies that enhance coherence between interacting units.

##  IX Acknowledgments

We thank Tiago Pereira for enlightening discussions and supporting to present
our results.

##  Appendix A Simulation of the expected time to the spark of synchronization

We present the algorithm we used for obtaining the expected time to
synchronization (Alg. 1). Its input consists of the adjacency matrix $A_{ij}$,
the network coupling strength $\alpha$, the states initial condition
$z^{0}_{i}$ (randomly sampled from a uniform distribution
${\mathcal{U}}_{[0,2\pi]}$), the maximum number of iterations $t_{max}$ and
the parameter $r_{c}$ that determines when the spark of synchronization
occurs. Because $r^{t}$ has different asymptotic behaviors for different
network configurations, to determined the value $r_{c}$ for each network
configuration, we first calculated the average value of $r$ after
synchronization and defined $r_{c}$ as half of it. We use half since the value
of $r$ at synchronization is not completely stable. For the simulations
presented in Figures 2 and 3, we considered $t_{max}=100,000$.

Algorithm 1 Time to emergence of synchrony

1:$z^{0}_{i},\alpha,A_{ij},r_{c},t_{max}$

2:$t_{s}$

3:$count\leftarrow 0$

4:for $t\in 0,...,t_{max}$ do

5: $\mbox{Compute }r^{t}$

6: if $r^{t}\geq r_{c}$ then

7: $count\leftarrow count+1$

8: $z^{t}\leftarrow N\mbox{ randomly sampled values from
}{\mathcal{U}}_{[0,2\pi]}$

9: else

10: $\mbox{Update }z^{t}$

11: end if

12:end for

13:$t_{s}\leftarrow t_{max}/count$

##  Appendix B Maps with connectivity $w$ close to the ratio $2/\alpha r^{t}$
are more likely to synchronize

Recall from Section II that
$V^{t}_{i}\sim{\mathcal{N}}^{2}(d_{i}V^{t},(d_{i}/2){\bf I})$. Thus, we can
approximate Equation 3 as follows,

$z^{t+1}_{i}\approx 2z^{t}_{i}+\alpha w_{i}\Im(V^{t}\bar{u}^{t}_{i}),$ (8)

with $V^{t}=r^{t}e^{i\theta^{t}}$. For all maps $i$, we can write
$z^{t}_{i}=\theta^{t}+h^{t}_{i}$, with $h^{t}_{i}$ small for those maps that
are synchronized at time $t$. Thus, by approximating Equation 8 by its first
degree (i.e., using $sin(h)\approx h$ for $h\approx 0$), we obtain
$z^{t+1}_{i}\approx 2\theta^{t}+h_{i}(2-w_{i}\alpha r^{t})$. The closest the
term $2-w_{i}\alpha r^{t}$ is to $0$, the more likely a map is to synchronize,
i.e., $z^{t+1}_{i}\approx 2\theta^{t}$. Therefore, maps with connectivity $w$
close to the ratio $2/\alpha r^{t}$ are more likely to synchronize.

##  Appendix C Proofs of the Lemmas

Recalling Definitions 3 and 5, here we present the proofs of Lemmas 4, 6, 7
and 8.

######  Lemma.

$E_{\tau_{a}}=a$.

###### Proof.

$\displaystyle E_{\tau_{a}}$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}e^{iz}\tau_{a}(z)dz$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{iz}(1+2\epsilon\cos(z-\theta))dz$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{iz}dz+\frac{\epsilon}{\pi}\int_{-\pi}^{\pi}e^{i(z^{\prime}+\theta)}\cos(z^{\prime})dz^{\prime}$
$\displaystyle=$
$\displaystyle\frac{\epsilon}{\pi}e^{i\theta}\int_{-\pi}^{\pi}e^{iz^{\prime}}\cos(z^{\prime})dz^{\prime}$
$\displaystyle=$
$\displaystyle\frac{\epsilon}{\pi}e^{i\theta}\left(\int_{-\pi}^{\pi}\cos(z^{\prime})^{2}dz^{\prime}+i\int_{-\pi}^{\pi}\sin(z^{\prime})\cos(z^{\prime})dz^{\prime}\right)$
$\displaystyle=$
$\displaystyle\frac{\epsilon}{\pi}e^{i\theta}\big{(}\pi+0\big{)}$
$\displaystyle=$ $\displaystyle\epsilon e^{i\theta}=a\qed$

######  Lemma.

Given $a\in{\mathbb{C}}$ and $\nu\in{\mathbb{R}}^{+}$,

$\mathbb{D}_{\nu}(\tau_{a})\approx\tau_{a^{\prime}},$

where $a^{\prime}=e^{-\nu^{2}/2}a$.

###### Proof.

Let $G_{\nu}(x)$ be the Gaussian density function
$\frac{1}{\nu\sqrt{2\pi}}e^{-\frac{1}{2\nu^{2}}x^{2}}$. We can write,

$\mathbb{D}_{\nu}(\tau)(z)=\int_{\mathbb{R}}\tau(z+\epsilon)G_{\nu}(\epsilon)d\epsilon.$

Then, using that
$\int_{\mathbb{R}}e^{i\epsilon}G_{\nu}(\epsilon)d\epsilon=e^{-\nu^{2}/2}$ and
writing $\tau_{a}(z)$ as $1+2\epsilon\Re(e^{i(z-\theta)})$, we obtain,

$\displaystyle\mathbb{D}_{\nu}(\tau_{a})(z)$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\tau_{a}(z+\epsilon)G_{\nu}(\epsilon)d\epsilon$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\left(1+2\epsilon\Re\left(e^{i(z-\theta+\epsilon)}\right)\right)G_{\nu}(\epsilon)d\epsilon$
$\displaystyle=$ $\displaystyle
1+2\epsilon\Re\left(\int_{\mathbb{R}}e^{i(z-\theta+\epsilon)}G_{\nu}(\epsilon)d\epsilon\right)$
$\displaystyle=$ $\displaystyle
1+2\epsilon\Re\left(e^{i(z-\theta)}\int_{\mathbb{R}}e^{i\epsilon}G_{\nu}(\epsilon)d\epsilon\right)$
$\displaystyle=$ $\displaystyle
1+2\epsilon\Re\left(e^{i(z-\theta)}\right)e^{-\nu^{2}/2}$ $\displaystyle=$
$\displaystyle\tau_{e^{-\nu^{2}/2}a}(z)\qed$

######  Lemma.

Given $a\in{\mathbb{C}}$ and $\nu\in{\mathbb{R}}^{+}$ small, $E_{\tau_{a}}=a$,
$E_{\mathbb{D}_{\nu}(\tau_{a})}=e^{-\nu^{2}/2}a$.

This lemma is immediate from the two previous lemmas.

######  Lemma.

Given $a\in{\mathbb{C}}$ and $P\in{\mathbb{C}}$ ,

${\mathbb{M}}_{P}(\tau_{a})\approx\tau_{a^{\prime}},$

where $a^{\prime}=P(P+4a)/8$.

###### Proof.

Let $a=(\epsilon/2)e^{i\delta}$, and $P=2me^{i\theta}$, where the factors of
$1/2$ and $2$ are used to simplify the notation below. By considering a
rotation, we may assume $\theta=0$ and $\delta=\delta-\theta$. This is evident
since, if $z^{\prime}=z-\theta$, then
$2z+m\sin(\theta-z)=2z^{\prime}+2\theta+m\sin(-z^{\prime})$. Thus, if we
define $\tau^{\prime}(z^{\prime})=\tau(z^{\prime}+\theta)$, then
${\mathbb{M}}_{P}(\tau)(z^{\prime}+2\theta)={\mathbb{M}}_{e^{-i\theta}P}(\tau^{\prime})(z^{\prime})$.
We will use this to find the formula for $\theta\neq 0$ later.

We are applying the function $f(z)=2z-2m\sin(z)$ to the distribution
$\tau(z)=(1/2\pi)(1+\epsilon\cos(z-\delta))$. Thus,
${\mathbb{M}}_{P}(\tau)(x)=\sum_{z\in
f^{-1}(x)}\frac{\tau(z)}{f^{\prime}(z)}$. For $m$ small, the function $f$ is
very close to $2z$, and hence, it is 2-to-1 as a function on $S^{1}$. If
$f(z)=x$, $z=x/2+m\sin(z)$, which we can approximate to the first order (using
$z\approx x/2$ and ignoring terms that contain $m^{2}$) by,

$z=x/2+m\sin(x/2)$

and

$2\pi\tau(z)=1+\epsilon\cos(x/2+m\sin(x/2)-\delta),$

which, using the 1st-order Taylor expansion of $cosine$ around $x/2-\delta$,
we can approximate to second order as,

$2\pi\tau(z)=1+\epsilon\cos(x/2-\delta)-\epsilon m\sin(x/2)\sin(x/2-\delta).$

The derivative of $f$ is given by $f^{\prime}(z)=2-2m\cos(z)$, which we can
approximate to second order by

$f^{\prime}(z)\approx 2-2m\cos(x/2+m\sin(x/2)).$

Considering again the first order Taylor expansion of $cosine$, we obtain,

$f^{\prime}(z)\approx 2-2m\cos(x/2)+2m^{2}\sin^{2}(x/2).$

Using that the second degree approximation to the inverse of $1+am+bm^{2}$ is
$1-am+(a^{2}-b)m^{2}$, we obtain,

$\displaystyle 2/f^{\prime}(z)$ $\displaystyle\approx$ $\displaystyle
1+m\cos(x/2)+m^{2}(\cos^{2}(x/2)-\sin^{2}(x/2))$ $\displaystyle\approx$
$\displaystyle 1+m\cos(x/2)+m^{2}\cos(x).$

There are two points in $f^{-1}(x)$, which we call $z_{1}$ and $z_{2}$. They
correspond to the fact that $x$ and $x+2\pi$ represent the same element in
$S^{1}$. So, we consider $z_{1}$ using $x$ as above and $z_{2}$ using
$x+2\pi$. Then, considering that $\cos(x/2+\pi)=-\cos(x/2)$ and
$\sin(x/2+\pi)=-\sin(x/2)$, we obtain,

$\displaystyle 2\pi\tau(z_{1})$ $\displaystyle=$ $\displaystyle 1-\epsilon
m\sin(x/2)\sin(x/2-\delta)+\epsilon\cos(x/2-\delta)$ $\displaystyle
2\pi\tau(z_{2})$ $\displaystyle=$ $\displaystyle 1-\epsilon
m\sin(x/2)\sin(x/2-\delta)-\epsilon\cos(x/2-\delta)$ $\displaystyle
2/f^{\prime}(z_{1})$ $\displaystyle\approx$ $\displaystyle
1+m^{2}\cos(x)+m\cos(x/2)$ $\displaystyle 2/f^{\prime}(z_{2})$
$\displaystyle\approx$ $\displaystyle 1+m^{2}\cos(x)-m\cos(x/2).$

Then, with the fact that $(A+B)(C+D)+(A-B)(C-D)=2AC+2BD$, we can obtain,

$\displaystyle 2\pi{\mathbb{M}}_{P}(\tau)(x)$ $\displaystyle=$ $\displaystyle
2\pi\tau(z_{1})/f^{\prime}(z_{1})+2\pi\tau(z_{2})/f^{\prime}(z_{2})$
$\displaystyle\approx$ $\displaystyle(1-\epsilon
m\sin(x/2)\sin(x/2-\delta))(1+m^{2}\cos(x))+\epsilon
m\cos(x/2-\delta)\cos(x/2)$ $\displaystyle\approx$ $\displaystyle 1-\epsilon
m\sin(x/2)\sin(x/2-\delta)+m^{2}\cos(x)+\epsilon m\cos(x/2-\delta)\cos(x/2)$
$\displaystyle=$ $\displaystyle 1+m^{2}\cos(x)+\epsilon
m\big{(}\cos(x/2-\delta)\cos(x/2)-\sin(x/2)\sin(x/2-\delta)\big{)}$
$\displaystyle=$ $\displaystyle 1+m^{2}\cos(x)+\epsilon m\cos(x-\delta)$
$\displaystyle=$ $\displaystyle 1+\Re(e^{ix}(m^{2}+\epsilon me^{-i\delta}))$
$\displaystyle=$ $\displaystyle 2\pi\tau_{b^{\prime}}(x),$

where $b^{\prime}=(m^{2}+\epsilon me^{-i\delta})/2$.

Now we remove the assumption that $\theta=0$. Recall that
${\mathbb{M}}_{P}(\tau)(z^{\prime}+2\theta)={\mathbb{M}}_{e^{-i\theta}P}(\tau^{\prime})(z^{\prime})$.
We thus need to move our result by $2\theta$, i.e.,

$\displaystyle 2\pi{\mathbb{M}}_{P}(\tau)(x)$ $\displaystyle=$ $\displaystyle
1+\Re(e^{i(x-2\theta)}(m^{2}+\epsilon me^{-i(\delta+\theta)}))$
$\displaystyle=$ $\displaystyle 1+\Re(e^{ix}\times
me^{-i\theta}\times(me^{-i\theta}+\epsilon e^{-i\delta}))$ $\displaystyle=$
$\displaystyle 2\pi\tau_{a^{\prime}}$

where $a^{\prime}=(P/2)(P/2+2a)/2=(P(P+4a))/8$. ∎

##  Appendix D The term $a_{w}$

In Section VI, the expression for the expected value of $V^{+}$ (Equation 7)
presents the term,

$\int w^{2}a_{w}e^{-\alpha^{2}w/4C}\delta(w),$ (9)

which contains the randomly distributed values $a_{w}$. Since $a_{w}$ is the
mean field of $\rho_{w}$, which is close to uniform before synchronization, we
assume that, for each $w$, $a_{w}$ is normally distributed around $0$ with
2D-variance $1/(2N\delta(w)dw)$. Notice that $N\delta(w)dw$ is the number of
maps with connectivity between $w$ and $w+dw$. Notice also that,

$\int wa_{w}\delta(w)dw=V.$ (10)

Now, for each $w$, consider a collection of a random variable $Z_{w}$ in the
complex plane, normally distributed with mean 0 and variance 1. Under these
circumstances, we can assume that $a_{w}$’s are described by
$Z_{w}/\sqrt{2N\delta(w)}$. Hereupon, consider $a$ and $Z$ as random variables
in ${\mathbb{R}}^{\mathbb{R}}$, with $a=T(Z)$, where $T$ is the map from
${\mathbb{R}}^{\mathbb{R}}\to{\mathbb{R}}^{\mathbb{R}}$ given by
$T(g)(w)=g(w)/\sqrt{2N\delta(w)}$. Let $C,D\in{\mathbb{R}}^{\mathbb{R}}$ be
the maps,

  * •

$C(w)=w\delta(w)$,

  * •

$D(w)=w^{2}e^{-w/4C}\delta(w).$

The motivation to introduce these maps is to write Equation 9 as,

$\langle D,a\rangle,$

where the inner product is defined as $\langle g,f\rangle=\int f(w)g(w)dw$.
This formulation also allow us to rewrite Equation 10 as,

$\langle C,a\rangle=V.$

Notice that $\langle f,T(g)\rangle=\langle T(f),g\rangle$ for all
$f,g\in{\mathbb{R}}^{\mathbb{R}}$. Therefore, we can procedure as follows: We
want to obtain the expected value of $\langle T(D),Z\rangle$, given that
$\langle T(C),Z\rangle=V$. To this end, we decompose $T(D)$ into orthogonal
functions as,

$T(D)=\beta T(C)+E$

where $T(C)$ and $E$ are orthogonal (i.e., $\langle T(C),E\rangle=0)$. By
doing so, we obtain,

$\langle D,a\rangle=\langle T(D),Z\rangle=\beta V+\langle E,Z\rangle.$

Thus, $\langle D,a\rangle$ has expected value $\beta V$, and variance
$||E||^{2}$. We calculate $\beta$ and $E$ as follows,

$\beta=\frac{\langle T(C),T(D)\rangle}{||T(C)||^{2}}\mbox{,}\quad E=T(D)-\beta
T(C).$

The numerator and denominator of $\beta$, respectively, are then obtained as,

$\displaystyle\langle T(C),T(D)\rangle$ $\displaystyle=$
$\displaystyle\int\frac{w\delta(w)}{\sqrt{2N\delta(w)}}\frac{w^{2}e^{-\alpha^{2}w/4C}\delta(w)}{\sqrt{2N\delta(w)}}dw$
$\displaystyle=$ $\displaystyle\frac{1}{2N}\int
w^{3}e^{-\alpha^{2}w/4C}\delta(w)dw$ $\displaystyle=$ $\displaystyle
K^{\alpha}_{\delta}(C),$ $\displaystyle\langle T(C),T(C)\rangle$
$\displaystyle=$
$\displaystyle\int\frac{w\delta(w)}{\sqrt{2N\delta(w)}}\frac{w\delta(w)}{\sqrt{2N\delta(w)}}dw$
$\displaystyle=$ $\displaystyle\frac{1}{2N}\int w^{2}\delta(w)dw$
$\displaystyle=$ $\displaystyle M_{2}.$

We conclude that $\int w^{2}a_{w}e^{-\alpha^{2}w/4C}\delta(w)$ has mean
$VK^{\alpha}_{\delta}(C)/M_{2}$, as needed in the derivation of Section VII.
As for $||E||^{2}$, we obtain,

$\langle
E,E\rangle=\frac{1}{2N}\int\big{(}w^{2}e^{-\alpha^{2}w/4C}-w(K^{\alpha}_{\delta}(C)/M_{2})\big{)}^{2}\delta(w)dw.$

This value is rather negligible in all examples we considered in this work.
For all 81 settings considered in Figure 2, the mean value for $\alpha||E||$
was 0.03 with maximum 0.07. Putting the results above with Equation 7, we
obtain that the randomness of $a_{w}$’s add 2D-noise with standard deviation
$||V||\alpha||E||$. Therefore, it is only a small percentage of the size of
$||V||$ that is being added as noise. The effect is probably not completely
negligible and more investigation is suggested for future works in this area.

## References

  * Milgram (1967) S. Milgram, Psychology today 2, 60 (1967).
  * Barabási and Albert (1999) A.-L. Barabási and R. Albert, science 286, 509 (1999).
  * Newman (2002) M. E. Newman, Physical review letters 89, 208701 (2002).
  * Radicchi et al. (2004) F. Radicchi, C. Castellano, F. Cecconi, V. Loreto, and D. Parisi, Proceedings of the national academy of sciences 101, 2658 (2004).
  * Zhou and Mondragón (2004) S. Zhou and R. J. Mondragón, IEEE communications letters 8, 180 (2004).
  * Barabasi and Oltvai (2004) A.-L. Barabasi and Z. N. Oltvai, Nature reviews genetics 5, 101 (2004).
  * Pastor-Satorras et al. (2001) R. Pastor-Satorras, A. Vázquez, and A. Vespignani, Physical review letters 87, 258701 (2001).
  * Arenas et al. (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports 469, 93 (2008), ISSN 0370-1573.
  * Sarfati et al. (2023) R. Sarfati, K. Joshi, O. Martin, J. C. Hayes, S. Iyer-Biswas, and O. Peleg, eLife 12, e78908 (2023), ISSN 2050-084X.
  * Winfree (2002) A. T. Winfree, Science 298, 2336 (2002), eprint https://www.science.org/doi/pdf/10.1126/science.1072560.
  * Strogatz (2004) S. Strogatz, _Sync: How Order Emerges From Chaos In the Universe, Nature, and Daily Life_ (Hyperion Books, 2004), ISBN 9780786887217.
  * Rodrigues et al. (2016) F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, Physics Reports 610, 1 (2016), ISSN 0370-1573, the Kuramoto model in complex networks.
  * Kuramoto (2013) Y. Kuramoto, _Chemical Oscillations, Waves, and Turbulence_ , Dover Books on Chemistry Series (Dover Publications, Incorporated, 2013), ISBN 9780486788593.
  * Lee (2005) D.-S. Lee, Phys. Rev. E 72, 026208 (2005).
  * Ichinomiya (2004) T. Ichinomiya, Phys. Rev. E 70, 026116 (2004).
  * Gómez-Gardeñes et al. (2007) J. Gómez-Gardeñes, Y. Moreno, and A. Arenas, Phys. Rev. E 75, 066106 (2007).
  * Baek and Ott (2004) S.-J. Baek and E. Ott, Phys. Rev. E 69, 066210 (2004).
  * Viana et al. (2005) R. Viana, C. Grebogi, S. de S.Pinto, S. Lopes, A. Batista, and J. Kurths, Physica D: Nonlinear Phenomena 206, 94 (2005), ISSN 0167-2789.
  * Pereira et al. (2017) T. Pereira, S. Strien, and M. Tanzi, Journal of the European Mathematical Society 22 (2017).
  * Corder et al. (2023) R. M. Corder, Z. Bian, T. Pereira, and A. Montalbán, Chaos: An Interdisciplinary Journal of Nonlinear Science 33, 091103 (2023), ISSN 1054-1500.
  * Komarov and Pikovsky (2015) M. Komarov and A. Pikovsky, Phys. Rev. E 92, 020901 (2015), URL https://link.aps.org/doi/10.1103/PhysRevE.92.020901.
  * Ottino-Löffler and Strogatz (2018) B. Ottino-Löffler and S. H. Strogatz, Phys. Rev. Lett. 120, 264102 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.264102.
  * Mendonca et al. (2023) H. M. Mendonca, R. Tönjes, and T. Pereira, Entropy 25 (2023), ISSN 1099-4300, URL https://www.mdpi.com/1099-4300/25/7/983.
  * Liang and Cvitanović (2022) H. Liang and P. Cvitanović, Journal of Physics A: Mathematical and Theoretical 55, 304002 (2022).
  * Stankovski et al. (2017) T. Stankovski, T. Pereira, P. V. E. McClintock, and A. Stefanovska, Rev. Mod. Phys. 89, 045001 (2017).
  * Restrepo et al. (2006) J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. Lett. 96, 254103 (2006).
  * Iaccarino et al. (2016) H. F. Iaccarino, A. C. Singer, A. J. Martorell, A. Rudenko, F. Gao, T. Z. Gillingham, H. Mathys, J. Seo, O. Kritskiy, F. Abdurrob, et al., Nature 540, 230 (2016).
  * Martorell et al. (2019) A. Martorell, A. Paulson, H.-J. Suk, F. Abdurrob, G. Drummond, W. Guan, J. Young, D. Kim, O. Kritskiy, S. Barker, et al., Cell 177 (2019).
  * Motter et al. (2013) A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Nature Physics 9, 191 (2013).