# Exponentially long transient time to synchronization of coupled chaotic
circle maps in dense random networks

Hans Muller Mendonca  Instituto de Ciências Matemáticas e Computação,
Universidade de São Paulo, São Carlos, São Paulo, Brazil  Ralf Tönjes
Institute of Physics and Astronomy, Potsdam University, 14476 Potsdam-Golm,
Germany  Tiago Pereira  Instituto de Ciências Matemáticas e Computação,
Universidade de São Paulo, São Carlos, São Paulo, Brazil

###### Abstract

We study the transition to synchronization in large, dense networks of chaotic
circle maps, where an exact solution of the mean-field dynamics in the
infinite network and all-to-all coupling limit is known. In dense networks of
finite size and link probability of smaller than one, the incoherent state is
meta-stable for coupling strengths that are larger than the mean-field
critical coupling. We observe chaotic transients with exponentially
distributed escape times and study the scaling behavior of the mean time to
synchronization.

synchronization; random networks; chaotic maps; mean-field analysis; finite
size effects

††preprint: APS/123-QED

##  I Introduction

Complex nonlinear systems often exhibit collective synchronization phenomena
which can play an important role for the overall functioning of a system
Lehnertz  _et al._ (2009); Schnitzler and Gross (2005); Pikovsky  _et al._
(2001). Phase oscillator models can elucidate key aspects of the mechanism
that generates the collective motion Stankovski  _et al._ (2017). The Kuramoto
model, for instance, is particularly useful in describing groups of weakly
coupled oscillators such as Josephson junctions, and they can be analyzed in
almost full detail in the thermodynamic limit of infinitely many oscillators.
Indeed, Kuramoto himself initially studied the fully connected networks of
coupled oscillators with frequency heterogeneity, and obtained the critical
value of the coupling strength for the transition from incoherence to
synchronized collective oscillations. Kuramoto (1984).

While such predictions are obtained in the thermodynamic limit, they have been
used as fruitful approaches to describe networks with finitely many
oscillators Bick  _et al._ (2020); Eldering  _et al._ (2021). However, recent
work has shown that finite size fluctuations or sparse connections in the
network can significantly impact on the overall dynamics. In fact, in certain
models, synchronization cannot, even approximately, be predicted from the
mean-field approximation in the thermodynamic limit Tönjes  _et al._ (2010).
That is, in these models, a transition to synchronization occurs or is
inhibited because of finite size fluctuations Komarov and Pikovsky (2015);
Ottino-Löffler and Strogatz (2018). The interplay between mean-field
predictions and finite-size fluctuations for general models remains elusive
and requires further investigation.

In this work, we study chaotic phase maps in dense networks where the mean-
field dynamics can be analyzed exactly in the thermodynamic limit. For small
coupling, due to the chaotic phase dynamics, only incoherence is stable. For a
range of coupling strengths, mean-field analysis predicts coexistence between
complete chaotic synchronization and incoherence, and for strong coupling, the
incoherence becomes unstable. Then, complete synchronization is the globally
attracting state in our model. Our results are two-fold:

(i) For coupling strengths with a stable coexistence of incoherence and
synchronization, although incoherence is locally attracting, finite-size
fluctuations can take the system into the basin of attraction of the absorbing
state of complete synchronization. Starting near incoherence with uniformly
distributed random oscillator phases, the distribution of transient times
towards synchronization is exponential and scales as a power of the system
size.

(ii) Above the critical coupling strength, in dense but incomplete networks,
although linear stability analysis of the mean-field equations suggests that
any nonzero mean field, e.g., finite size fluctuations of the mean field, will
grow exponentially fast, we observe an exponentially long chaotic transient in
the incoherent state. Such a delayed transition to synchronization has so far
not been described in dense networks of coupled phase oscillators or coupled
chaotic maps.

##  II Model of Coupled Chaotic Maps

The local phase dynamics in each node is modelled as a Bernoulli map of the
circle with time steps $t\in\mathbb{Z}$

$\varphi(t+1)=f(\varphi(t))=2\varphi(t)\mod 2\pi,$ (1)

or via the abuse of notation on the complex unit circle $z=\exp(i\varphi)$, we
write $z(t+1)=f(z(t))=z(t)^{2}$. This map is chaotic and structurally stable
Tanzi  _et al._ (2019). That is, the statistical properties of the map persist
under small perturbations. Therefore, for small coupling, the maps behave as
nearly independent, and no collective dynamics is possible for small coupling.
In Gong  _et al._ (2020), the global coupling of the phase dynamics is
implemented as a Moebius map on the complex unit circle. The Moebius map has
been shown to give exact solutions of sinusoidally forced phase dynamics
Marvel  _et al._ (2009), including the Kuramoto model, Winfree-type phase
equations, and via a nonlinear transformation, the dynamics of theta neurons
Montbrió  _et al._ (2015). It is therefore a meaningful alternative to the
sine coupling in the standard circle map. Here, we use a composition of (1)
and a Moebius map (see Figure 1)

$z{(t+1)}=M\left(f(z(t)),\Phi(t),\tau(t)\right),$ (2)

where

$M(w,\Phi,\tau)=\frac{e^{i\Phi}\tau+w}{1+e^{-i\Phi}\tau w}$ (3)

for a coupling intensity $-1<\tau<1$, an angle of contraction
$\Phi\in\mathbb{S}^{1}$, and a point
$w\in\mathbb{D}=\\{z\in\mathbb{C}:|z|<1\\}$ on the open complex unit disc. The
family of Moebius maps is a group of biholomorphic automorphisms of
$\mathbb{D}$, and via analytic continuation, these transformations map the
boundary of $\mathbb{D}$ bijectively onto itself. The effect of (3) on the
unit disc is a contraction of almost all points towards $\exp(i\Phi)$ on the
boundary where $\lim_{\tau\to\pm 1}M(w,\Phi,\tau)=\pm\exp(i\Phi)$ and
$\lim_{\tau\to 0}M(w,\Phi,\tau)=w$. The parameter $\tau$ characterizes the
strength of the contraction. For $\tau\to 0$, the map (2) approaches the
uncoupled dynamics (1). Moreover, the family of wrapped Cauchy distributions

$p(\varphi)=\frac{1}{2\pi}\frac{1-R^{2}}{|1-Re^{i(\varphi-\Theta)}|^{2}}$ (4)

which includes incoherence as the uniform distribution when $R\to 0$ and a
delta distribution at $\varphi=\Theta$ when $R\to 1$, is invariant under (2)
and (3) Gong  _et al._ (2020); Marvel  _et al._ (2009); Pikovsky (2021). This
family of continuous phase measures, in the context of phase synchronization,
is known as the Ott-Antonsen manifold, and assuming this form of phase
distribution is equivalent to the so called Ott-Antonsen ansatz Ott and
Antonsen (2008, 2009). The Ott-Antonsen manifold is parameterized using the
mean-field amplitude $R$ and the mean-field angle $\Theta$

$Z=Re^{i\Theta}=\int_{0}^{2\pi}e^{i\varphi}p(\varphi)\,d\varphi.$ (5)

The mean-field amplitude $R$ is the Kuramoto order parameter Kuramoto (1975),
which is zero for incoherence, i.e., a uniform phase distribution, and $R=1$
for complete synchronization $\varphi_{n}=\Theta$ (a.s.). Furthermore, the
higher circular moments $Z_{q}$ on the Ott-Antonsen manifold with
$q\in\mathbb{Z}$ are integer powers of the mean field

$\displaystyle Z_{q}=\int_{0}^{2\pi}e^{iq\varphi}p(\varphi)\,d\varphi=Z^{q}.$
(6)

As a consequence, phase doubling maps the circular moments as
$f(Z_{q}(t))=Z_{2q}(t)=Z_{2}^{q}(t)=f(Z_{1}(t))^{q}$, leaving the Ott-Antonsen
manifold invariant and mapping the mean-field amplitude and phase as $R\to
R^{2}$ and $\Theta\to 2\Theta$.

Figure 1: Dynamics of phases. $N=30$ points on the complex unit circle colored
by phase, and corresponding mean field (red dot inside the unit circle). From
left to right : initial phase configuration at the points $z_{n}$ with mean-
field amplitude $R=0.5$ and mean-field phase $\Theta=\pi/4$, after chaotic
phase doubling $z_{n}^{2}$ with $R^{2}=0.25$ and $2\Theta=\pi/2$, and after
subsequent contraction toward the angle $\pi/2$ with intensity $\tau=0.5$.

To couple the dynamics of the Bernoulli maps (2), the parameters $\Phi(t)$ and
$\tau(t)$ in (3) should be defined as functions of the ensemble mean field.
Following Gong  _et al._ (2020), we define the contraction angle $\Phi(t)$ and
the coupling intensity $\tau(t)$ as

$\displaystyle Z(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}z_{n}(t)=R(t)e^{i\Theta(t)}$ (7)
$\displaystyle\Phi(t)$ $\displaystyle=$ $\displaystyle 2\Theta(t)$ (8)
$\displaystyle\tau(t)$ $\displaystyle=$
$\displaystyle\tanh\left(\frac{\varepsilon}{2}R(t)\right),$ (9)

where $\varepsilon$ is a coupling strength. For $\tau=1$, when $\varepsilon
R\to\infty$, the phases are contracted to a single point $\exp(2i\Theta)$ on
the unit circle. For small values of $\varepsilon R$, we can expand (2) to the
linear order and obtain the more familiar form of mean-field coupled circle
maps with phase doubling

$\varphi_{n}(t+1)=2\varphi_{n}(t)+\varepsilon
R(t)\sin\left(2\Theta(t)-2\varphi_{n}(t)\right)+O(\varepsilon^{2}R^{2}(t)).$
(10)

The crucial observation is that on the Ott-Antonsen manifold, the mean-field
$Z=R\exp(i\Theta)$ transforms exactly the same way via (2),(3) as each element
$z=\exp(i\varphi)$ on the unit circle Marvel  _et al._ (2009); Gong  _et al._
(2020); that is,

$Z(t+1)=M(Z^{2}(t),\Phi(t),\tau(t)).$ (11)

It is highly unusual that a closed analytic expression for the dynamics of the
mean field can be derived and thus analyzed in coupled nonlinear dynamical
systems. The reduction in infinitely dimensional microscopic dynamics to the
low-dimensional dynamics of the mean-field Ott and Antonsen (2008) has been
tremendously successful in the analysis of synchronization phenomena over the
last decade, while the effects of the finite system size $N$ remain difficult
to analyze Peter and Pikovsky (2018); Peter  _et al._ (2019). We note that the
point measure of a finite ensemble of phases is never actually on the Ott-
Antonsen manifold, but can, in some sense, be arbitrarily close to the so-
called thermodynamic limit, i.e., the limit of the infinite system size
$N\to\infty$.

Applying the Ott-Antonsen ansatz to networks of phase oscillators is possible
if the network structure allows for the partitioning of the vertices into a
few classes of equivalent vertices. Assuming that all vertices of a class are
subjected to the same sinusoidal forcing, the dynamics of the phases in the
network can be reduced to the dynamics of coupled mean fields on the Ott-
Antonsen manifold for each vertex class Skardal and Restrepo (2012); Martens
(2010); Ottino-Löffler and Strogatz (2018); Martens  _et al._ (2016); Tönjes
_et al._ (2021). Additionally, heterogeneity in the oscillators and
fluctuations in the forces can be incorporated into the mean field dynamics if
they follow Cauchy distributions Laing (2009); Tönjes and Pikovsky (2020);
Clusella and Montbrió (2022).

###  II.1 Mean-Field Analysis

The mean-field dynamics (11) can be written in terms of the polar
representation

$\Theta(t+1)=f(\Theta(t))\,\,\,\,\,\mbox{~{}and
}\,\,\,\,\,R(t+1)=\frac{\tau(t)+R^{2}(t)}{1+\tau(t)R^{2}(t)}.$ (12)

This means that the dynamics of the phase $\Theta$ decouples from the
amplitude and will evolve chaotically. Using Equations (9) and (12), we obtain
the amplitude dynamics

$R(t+1)=\frac{\tanh\left(\frac{1}{2}\varepsilon
R(t)\right)+R^{2}(t)}{1+\tanh\left(\frac{1}{2}\varepsilon
R(t)\right)R^{2}(t)}$ (13)

which describes the exact evolution of the order parameter $R$ in a closed
form. We can readily determine the fixed points of the mean-field amplitude
$R(t)$ and their linear stability. Both the complete synchronization $R=1$ and
the complete desynchronization $R=0$ are fixed points of (13), and change
stability at unique critical points $\varepsilon_{1}=\ln(2)\approx 0.69$ and
$\varepsilon_{0}=2$, respectively, as determined by the eigenvalues of
Jacobian of Equation (13) at these fixed points. These critical points are
connected by an unstable fixed point branch $(\varepsilon(R_{u}),R_{u})$,
where

$\varepsilon(R_{u})=\frac{1}{R_{u}}\log\left(\frac{(1+R_{u})^{2}}{1+R_{u}^{2}}\right).$
(14)

This expression is derived from (13) by setting $R(t+1)=R(t)=R_{u}$ and
resolving the equation for $\varepsilon$.

This means that this system of all-to-all coupled, identical chaotic phase
maps will always evolve to complete synchronization or complete
desynchronization, with a small region $\ln(2)<\varepsilon<2$ of bistability
(Figure 2a).

###  II.2 Extension to Networks

Next, we have studied the same phase dynamics on a random network of $N$ maps
which are coupled to exactly $k$ different, random neighbors. Here, each phase
$\varphi_{n}$ couples to a local mean field

$Q_{n}=R_{n}e^{i\Theta_{n}}=\frac{1}{k}\sum_{n=1}^{N}A_{nm}z_{m}$ (15)

where $A_{nm}$ are the entries of the adjacency matrix, i.e., equal to one if
there is a link from vertex $m$ to vertex $n$, but zero otherwise, and
$k=\sum_{m=1}^{N}A_{nm}$ is the in-degree of node $n$, which, for
computational simplicity, we assume to be identical for all nodes. Thus, with
$\tau_{n}=\tanh\left(\frac{\varepsilon}{2}R_{n}\right)$, the dynamics of the
phases coupled through a network are

$z_{n}(t+1)=\frac{e^{2i\Theta_{n}(t)}\tau_{n}(t)+z_{n}^{2}(t)}{1+e^{-2i\Theta_{n}(t)}\tau_{n}(t)z_{n}^{2}(t)}.$
(16)

A class of networks is dense if $\lim_{N\to\infty}\langle k\rangle/N=p>0$,
where $\langle k\rangle$ is the mean node degree. Therefore, $p$ is the
fraction of nodes, in relation to the system size $N$, that an oscillator is
coupled to. Since dense networks are defined in the limit of $N\to\infty$,
there is no sharp distinction between sparse and dense networks of finite
size. We refer to a finite network as dense if two nodes share more than one
neighbor on average, i.e., $\langle k\rangle^{2}/N=p^{2}N>1$. In large dense
networks, the local mean fields of the oscillators in the neighborhood of each
node (15) are equal to the global mean field, with a deviation of
$O(1/\sqrt{k})$, where $k$ is the size of the neighborhood, i.e., the in-
degree of the node. Therefore, mean-field theory should be exact for dense
networks in the thermodynamic limit where $\langle k\rangle\to\infty$.

The network model First, we wish to compare the simulation results directly
with our mean-field analysis. For large random networks with a link density
$p=k/N$ and $0<p<1$, the numerical simulations are time-consuming since the
$N$ local mean fields at each node in the network need to be computed in each
time step. To simplify these computations, we use a random network where each
node couples to exactly $k$ different random neighbors. This model with a
unique in-degree of $k$ for each node is slightly different from the Erdös
Renyi model, with a Poissonian in-degree distribution of small relative width
$\textrm{std}(k)/\langle k\rangle\sim 1/\sqrt{k}$. For large $k$, the results
of the simulations in our random network model and other random networks with
uncorrelated node degrees and a vanishing relative width of the degree
distribution are expected to be identical.

(a)

(b)

Figure 2: Bifurcation diagram of the mean-field amplitude and a representation
of the network interaction. In (a), the bifurcation diagram of the all-to-all
coupling mean-field dynamics (12), i.e., on the Ott-Antonsen manifold. Dotted
lines show linearly unstable fixed points and solid lines show linearly stable
fixed points in the thermodynamic limit. (b) Venn diagram of a dense network
with $N$ vertices and connection probability $p$. The sets of neighbors of
nodes $m$ and $n$ are of size $pN$ and their overlap is of size $p^{2}N$,
resulting in correlated local mean fields $Q_{m}=R_{m}\exp(i\Theta_{m})$ and
$Q_{n}=R_{n}\exp(i\Theta_{n})$ acting on the states $z_{m}$ and $z_{n}$. The
ratio of the amplitudes of the local mean fields and the global mean field are
independent of the network size $N$.

##  III Results

###  III.1 Distributions of Transient Times

We perform a large number $M$ of simulations $m=1\ldots M$ from independent,
uniformly distributed random initial phases over a maximum of $T$ steps and
record in each simulation the first time step $t_{m}$ when $R\geq 0.5$, i.e.,
the transition time from an incoherent state to complete synchronization.
Finite-size scaling for such a discontinuous transition is challenging de
Oliveira  _et al._ (2015). The exponential distribution of the times $t_{m}$,
according to some characteristic transition rate, can be checked in a rank
plot of time points $t_{m}$, which gives the sample complementary cumulative
distribution $C(t)=\textrm{prob}(t\geq t_{m})=\textrm{rank}(t_{m})/M$
LABEL:Fig:Fig04a,d).

An exponential tail distribution $C(t)$ up to observation time $T$ indicates
an exponential distribution of transient times. Since the simulation time is
finite, transition times $t_{m}\geq T$ are not observed, which represents a
problem when we are interested in the average time to synchronization.
However, assuming a discrete exponential, i.e., geometric distribution, a
maximum likelihood estimation of the average transition time is possible up to
values considerably exceeding the observation time $T$ (see Appendix
LABEL:app1).

Denoting the number of simulations that synchronize at times $t_{m}<T$ as
$M_{T}$, and defining the observable values $l_{m}=\textrm{min}(t_{m},T)$, the
maximum likelihood estimation of the expected value
$T_{esc}=\textrm{E}[t_{m}]$ for the geometric distribution is

$T_{esc}=\frac{\left\langle l_{m}\right\rangle M}{M_{T}}.$ (17)

with the sample mean $\langle l_{m}\rangle$. If the transition to
synchronization is observed in all simulations, i.e., $M_{T}=M$, the estimator
is simply the sample mean of $t_{m}$, which is an estimator of $T_{esc}$ for
arbitrary transient time distributions. However, when most runs do not
synchronize within the finite simulation time $T$, the ratio $M/M_{T}$
contains additional information, and the estimated mean escape time can be
much larger than the observation time.

###  III.2 During Coexistence: Escape over the Unstable Branch

In Luccioli  _et al._ (2012), it was reported that the transition from
incoherence to collective dynamics in sparse networks of coupled logistic maps
is of the mean-field type. The analysis in Restrepo  _et al._ (2005) predicts
a shift in the critical coupling strength in random networks of Kuramoto phase
oscillators of the order $\langle k\rangle^{2}/\langle k^{2}\rangle$ due to
degree inhomogeneity, and $1/\langle k\rangle$ due to finite size fluctuations
of the local mean fields. That is, in dense, homogeneous networks with
$\langle k\rangle^{2}/\langle k^{2}\rangle\to 1$ and $\langle
k\rangle\to\infty$, the critical coupling strength does not change. We
expected to find similar behaviors for network-coupled Bernoulli maps. In
complete or almost complete networks $k/N=p\approx 1$ for $\varepsilon<2$,
there is a small probability that finite size fluctuations bring the order
parameter $R$ above the unstable branch, leading to a spontaneous transition
to complete synchronization, as shown in Figure 3a. We first observe the
scaling of the transient time in fully connected networks with $p=1$. For
values of $\varepsilon<\varepsilon_{0}=2.0$, the transition rate to
synchronization scales strongly with the size $N$ of the system (Figure
LABEL:Fig:Fig04b,c). However, for values $\varepsilon>\varepsilon_{0}$, the
average transition time depends very weakly on $N$, as the system grows
exponentially fast from a state of incoherence, with $R\approx 1/\sqrt{N}$. We
estimate a finite size scaling exponent $\beta$ below the transition threshold
by collapsing the curves $T_{esc}(\varepsilon,N)$ using the ansatz
$T_{esc}(\varepsilon,N)=T_{esc}\left((\varepsilon-\varepsilon_{0})N^{\beta}\right).$
The data are consistent with an ad hoc exponent of $\beta=1/3$ (Figure
LABEL:Fig:Fig04c).

###  III.3 Above the Critical Coupling Strength: Long Chaotic Transient

Above the critical coupling strength $\varepsilon>\varepsilon_{0}=2$, we
expected finite size fluctuations to grow exponentially fast and independently
of $N$, as predicted by linear stability analysis of the mean-field equations
(13). Instead, for small connection probabilities $0<p<1$, we have observed a
chaotic transient with seemingly stationary finite size fluctuations
$O(1/\sqrt{N})$ of the mean field (Figure 3). In the large $N$ limit, the
distribution of the transient times depends on the link density $p$ with
increasingly long transients as $p$ is decreased, but it is otherwise
independent of $N$.

A coupling strength for which a transition to complete synchronization could
still be observed within the simulation time was considerably larger than the
mean-field critical coupling $\varepsilon_{0}=2$. That is, even in dense
networks and above the mean-field critical coupling, finite size fluctuations
will not necessarily result in the nucleation and exponential growth of a
collective mode. Such a delayed transition to synchronization Baer  _et al._
(1989) has so far not been described in systems of coupled phase oscillators
Ichinomiya (2004); Restrepo  _et al._ (2005); Ko and Ermentrout (2008) or
coupled logistic maps Luccioli  _et al._ (2012).

In Figure LABEL:Fig:Fig04f, we plot $T_{esc}$ over
$(\varepsilon-\varepsilon_{0})p$ to demonstrate that the average transition
time is roughly scaling as $1/p$. We do not look for higher-order corrections
such as a weak dependence of $\varepsilon_{0}$ on $p$, although the curves do
not collapse perfectly. Note that the escape time is largely independent of
the network size (Figure LABEL:Fig:Fig04e,f). For $p=0.1$, $0.05$, and $0.025$
we have performed simulations with $N=10^{4}$ (circles) and with $N=5\times
10^{4}$ (crosses) for comparison. For $p=0.01$, we compare network sizes
$N=10^{4}$ (circles) with very time-consuming simulations in networks with
$N=10^{5}$ (crosses).

(a)

(c)

(b)

(d)

Figure 3: Transient to synchronization for $N=$ 10,000 coupled maps in (a,b),
a fully connected network with coupling strength $\varepsilon=1.81$ below the
critical coupling $\varepsilon_{0}=2$, and (c,d), in a random network with
connection probability $p=0.1$ for a coupling strength of $\varepsilon=2.3$
above the critical coupling. The upper panels (a,c) show the order parameter
$R(t)$, and the lower panels (b,d), the real part of the ratio of the first
two circular moments $\textrm{Re}[Z_{1}^{2}/Z_{2}]$. This serves as a visual
measure of the alignment of the system state with the Ott-Antonsen manifold,
where the ratio is exactly equal to one. The dashed line in (a) marks the
value of the unstable fixed point of the mean-field dynamics, $R_{u}=0.098$.
Above that value, the state of complete synchronization is attractive on the
Ott-Antonsen manifold. In (b,d), the incoherent state $R=0$ is unstable;
however, finite size fluctuations do not grow exponentially. Instead, we
observe a long chaotic transient.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Statistics of transient times $t_{m}$ to synchronization. (a–c) In
the fully connected network; (d–f) in random networks of various link
densities $p=k/N$. The left panels show straight lines in semi-logarithmic
plots of cumulative tail distributions of the transient times, demonstrating
the rate character of the transition process. The middle panels show the
estimated average transient times for various combinations of system sizes
$N$, coupling strengths $\varepsilon$, and link densities $p$. The mean field
critical coupling strength $\varepsilon_{0}=2.0$ and the maximum observation
time $T$ are marked by dashed lines. In the globally coupled system in pannels
(a–c), the transient time depends strongly on the system size $N$, whereas in
dense networks and above $\varepsilon_{0}$ (d–f), the transient time depends
strongly on the link density $p=k/N$, but not on the system size. We
demonstrate the scaling of the transient times in panels (c) and (f) on the
right. In the globally coupled system, the exponential divergence of the
transient times below $\varepsilon_{0}$ appears to be a function of
$(\varepsilon-\varepsilon_{0})N^{\frac{1}{3}}$. In dense networks, the
exponential divergence is roughly a function of
$(\varepsilon-\varepsilon_{0})p$.

###  III.4 Discussion of Finite Size Scaling

Mean field theory assumes a phase distribution on the Ott-Antonsen manifold.
The characteristic function of a wrapped Cauchy distribution is the geometric
sequence $Z_{q}=Z^{q}$ of circular moments (6). However, in the incoherent
state with $N$ independent uniformly distributed phases $\varphi_{n}$ the
circular moments of an ensemble

$Z_{q}=\frac{1}{N}\sum_{n=1}^{N}e^{iq\varphi_{n}}$ (18)

are almost independent complex numbers with a Gaussian distribution of mean
zero and a variance of $1/N$ by virtue of the central limit theorem. The
action of the Bernoulli map on the circular moments is the shift

$Z_{q}\to Z_{2q},$ (19)

that is, it is achieved by discarding all odd circular moments. The
exponential growth of the order parameter in accordance to mean field theory
is expected after the distribution comes close to the Ott-Antonsen manifold,
i.e., when the first few circular moments align by chance sufficiently under
the mapping (19); in particular, $Z_{2}(t)\approx Z_{1}^{2}(t)$. Unless the
directions of $Z_{2}$ and $Z_{1}^{2}$ align by chance, as they would on the
Ott-Antonsen manifold, the subsequent contraction of strength $\varepsilon R$
in the direction of $Z_{1}^{2}$ after the phase doubling may even decrease the
amplitude of the order parameter. In addition, for coupling strengths
$\varepsilon$ below the critical value, $R=|Z_{1}|$ must be above the unstable
branch $R>R_{u}(\varepsilon)\sim(\varepsilon_{0}-\varepsilon)$.

The rate of such a random event should depend on the ratio between
$R_{u}(\varepsilon)$ and the standard deviation $1/\sqrt{N}$ of the Gaussian
distribution of the complex mean field. Based on this scaling argument, the
expected time to synchronize should scale
as$T_{esc}=T_{esc}((\varepsilon-\varepsilon_{0})\sqrt{N})$ below the critical
coupling. The best collapse of the estimated escape times in fully connected
networks of coupled Bernoulli maps was observed by scaling the distance to
$\varepsilon_{0}$ with $N^{1/3}$ (Figure 3c), i.e., the exponential divergence
of the escape time approaches $\varepsilon_{0}$ slower than $1/\sqrt{N}$ in
the thermodynamic limit. One possibility for this discrepancy is that the
scaling argument only considers the chance of $R>R_{u}$ and not the alignment
process of the higher-order circular moments.

Above the critical coupling strength, there is only the condition of the
alignment of circular moments with the Ott-Antonsen manifold for the
initiation of exponential growth. Since in the incoherent state, all circular
moments are random Gaussian with identical variance, the alignment process
(19) is strictly independent of the system size $N$. Once exponential growth
in the direction of the Ott-Antonsen manifold occurs, the time to
synchronization is logarithmic, that is, it is weakly dependent on $N$.
However, it appears that the alignment with the Ott-Antonsen manifold needs to
be stronger for networks with link densities of $p<1$. For small link
densities, the divergence of the escape time occurs at larger values
$\varepsilon>\varepsilon_{0}$. This is reminiscent of stabilization by noise
Khasminskii (2011), where a system is driven away from a low-dimensional
unstable manifold of a fixed point into stronger attracting stable directions.

In simulations of dense random networks of coupled Bernoulli maps, we could
see the independence of the mean escape time from the network size and the
scaling of the escape time with roughly $\sim 1/p$ (Figure LABEL:Fig:Fig04f).
To explain this scaling, we argue that mean field theory might be extended to
dense networks, where each node couples to a finite neighborhood of $pN$ nodes
in the network, and for every two nodes, these neighborhoods overlap on a set
of size $p^{2}N$ (Figure 2b). The local mean fields are Gaussian random forces
of mean value $Z$, variance $1/k=1/pN$, and a pairwise correlation of $p$,
which is the relative size of the overlap. The decrease in correlation between
the local mean fields in networks with link densities $p<1$ can be interpreted
as individual, finite size noise on the maps, which couple to the global mean
field, plus some uncorrelated random deviation. Therefore, the contractions of
the phases do not occur in the same direction for different nodes in the
network. The strength of the contraction in the direction of the mean field is
effectively reduced by the factor $p$, i.e.,

$\displaystyle\tau=\tanh\left(\frac{1}{2}\varepsilon
R\right)p\approx\frac{1}{2}\varepsilon pR$ (20)

shifting the coupling strength dependence of the transition time (above
$\varepsilon_{0}$) by a factor of $1/p$.

##  IV Conclusions

We have investigated the synchronization of coupled chaotic maps in dense
random networks, utilizing mean-field equations and examining network
configurations with different link probabilities. Firstly, we noticed the
existence of chaotic transients to synchronization within these networks. This
means that the incoherent state can persist for extended periods before
transitioning into synchronization. This finding led us to study the
statistics of transient times and their scaling behaviors in the process of
synchronization. The transition times follow exponential distributions,
indicating spontaneous transitions at a constant rate. It is noteworthy that
the transition from incoherence to complete synchronization only occurs
spontaneously in networks of finite size. Additionally, we have observed a
remarkable dependence of the transient times to synchronization on the link
probability $p$, represented by the ratio of the in-degree to the total number
of nodes, at coupling strengths where an immediate transition to synchrony
would be expected from mean field theory. Whether such a delayed transition is
due to the specifics of our model or is typical for a more general class of
dynamics remains an open question.

This research was funded by the FAPESP CEMEAI 391, Grant No. 2013/07375-0,
Serrapilheira Institute (Grant No.Serra-392 1709-16124), Newton Advanced
Fellow of the Royal Society (393 NAF$\backslash$R1$\backslash$180236), CAPES
and CNPq, Grant No 166191/2018-3.

##  V Apendix

Here, we calculate the maximum likelihood estimation for the mean value of a
geometric distribution $P(t;\alpha)=(1-\alpha)\alpha^{t}$ for discrete values
$t=0,1,\ldots$ of time steps when only times $t<T$ can be observed. The
expected value for the geometric distribution is

$\textrm{E}\left[t\right]=(1-\alpha)\sum_{t=0}^{\infty}t\alpha^{t}=\frac{\alpha}{1-\alpha}.$
(21)

Since the times $t_{m}$, $m=1\ldots M$ are only observable up to step $T-1$,
we define $l_{m}=\textrm{min}(t_{m},T)$. The probabilities for the possible
values of $l_{m}$ are

$\displaystyle P(l_{m}=T;\alpha)$ $\displaystyle=$ $\displaystyle
1-(1-\alpha)\sum_{t=0}^{T-1}\alpha^{t}=\alpha^{T}$ (22) $\displaystyle
P(l_{m}=t<T;\alpha)$ $\displaystyle=$ $\displaystyle(1-\alpha)\alpha^{t}.$
(23)

The derivative of the log-likelihood of $M$ independent observations $l_{m}$
with respect to the parameter $\alpha$ is

$\frac{\partial_{\alpha}P(l_{1},l_{2},\ldots,l_{M};\alpha)}{P(l_{1},l_{2},\ldots,l_{M};\alpha)}=\sum_{m=1}^{M}\frac{\partial_{\alpha}P(l_{m};\alpha)}{P(l_{m};\alpha)}.$
(24)

For the probabilities (22),(23), the derivatives are

$\displaystyle\frac{\partial_{\alpha}P(l_{m}=T,\alpha)}{P(l_{m}=T,\alpha)}$
$\displaystyle=$ $\displaystyle\frac{T}{\alpha}$ (25)
$\displaystyle\frac{\partial_{\alpha}P(l_{m}=t<T,\alpha)}{P(l_{m}=t<T,\alpha)}$
$\displaystyle=$ $\displaystyle\frac{t}{\alpha}-\frac{1}{1-\alpha}.$ (26)

For a maximum of the log-likelihood for the observed values $l_{m}$, the sum
in (24) is required to be zero. Inserting $M-M_{T}$ times the term (25) for
all observations $l_{m}=T$ and $M_{T}$ terms (26), one for each observation
$l_{m}=t<T$, we obtain

$(M-M_{T})\frac{T}{\alpha}+\sum_{l_{m}<T}\frac{l_{m}}{\alpha}-M_{T}\frac{1}{1-\alpha}=0.$
(27)

With

$\langle
l_{m}\rangle=\frac{1}{M}\sum_{m=1}^{M}l_{m}=\frac{1}{M}\left((M-M_{T})T+\sum_{l_{m}<T}l_{m}\right)$
(28)

we can divide (27) by the number $M$ of observations and re-order the equation
to obtain

$\frac{\langle l_{m}\rangle M}{M_{T}}=\frac{\alpha}{1-\alpha}.$ (29)

However, this is exactly the expected value $\textrm{E}\left[t\right]$ of time
steps for the full geometric distribution (21).

## References

  * Lehnertz  _et al._ (2009) K. Lehnertz, S. Bialonski, M.-T. Horstmann, D. Krug, A. Rothkegel, M. Staniek, and T. Wagner, Journal of neuroscience methods 183, 42 (2009).
  * Schnitzler and Gross (2005) A. Schnitzler and J. Gross, Nature reviews neuroscience 6, 285 (2005).
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, and J. Kurths, Self 2, 3 (2001).
  * Stankovski  _et al._ (2017) T. Stankovski, T. Pereira, P. V. E. McClintock, and A. Stefanovska, Rev. Mod. Phys. 89, 045001 (2017).
  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves and Turbulence_, Vol. 19 (Springer-Verlag, 1984).
  * Bick  _et al._ (2020) C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, The Journal of Mathematical Neuroscience 10, 9 (2020).
  * Eldering  _et al._ (2021) J. Eldering, J. S. Lamb, T. Pereira, and E. R. dos Santos, Nonlinearity 34, 5344 (2021).
  * Tönjes  _et al._ (2010) R. Tönjes, N. Masuda, and H. Kori, Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 033108 (2010).
  * Komarov and Pikovsky (2015) M. Komarov and A. Pikovsky, Phys. Rev. E 92, 020901 (2015).
  * Ottino-Löffler and Strogatz (2018) B. Ottino-Löffler and S. H. Strogatz, Phys. Rev. Lett. 120, 264102 (2018).
  * Tanzi  _et al._ (2019) M. Tanzi, T. Pereira, and S. van Strien, Ergodic Theory and Dynamical Systems 39, 1121 (2019).
  * Gong  _et al._ (2020) C. C. Gong, R. Toenjes, and A. Pikovsky, Physical Review E 102, 022206 (2020).
  * Marvel  _et al._ (2009) S. A. Marvel, R. E. Mirollo, and S. H. Strogatz, Chaos 19, 043104 (2009), arXiv:0904.1680 [nlin.SI] .
  * Montbrió  _et al._ (2015) E. Montbrió, D. Pazó, and A. Roxin, Phys. Rev. X 5, 021028 (2015).
  * Pikovsky (2021) A. S. Pikovsky, Applied Nonlinear Dynamics 29, 78 (2021).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 037113 (2008).
  * Ott and Antonsen (2009) E. Ott and T. M. Antonsen, Chaos: An interdisciplinary journal of nonlinear science 19, 023117 (2009).
  * Kuramoto (1975) Y. Kuramoto, in  _International symposium on mathematical problems in theoretical physics_ (Springer, 1975) pp. 420–422.
  * Peter and Pikovsky (2018) F. Peter and A. Pikovsky, Phys. Rev. E 97, 032310 (2018).
  * Peter  _et al._ (2019) F. Peter, C. C. Gong, and A. Pikovsky, Phys. Rev. E 100, 032210 (2019).
  * Skardal and Restrepo (2012) P. S. Skardal and J. G. Restrepo, Physical Review E 85, 016208 (2012).
  * Martens (2010) E. A. Martens, Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 043122 (2010).
  * Martens  _et al._ (2016) E. A. Martens, M. J. Panaggio, and D. M. Abrams, New Journal of Physics 18, 022002 (2016).
  * Tönjes  _et al._ (2021) R. Tönjes, C. E. Fiore, and T. Pereira, Nature Communications 12, 72 (2021).
  * Laing (2009) C. R. Laing, Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 013113 (2009).
  * Tönjes and Pikovsky (2020) R. Tönjes and A. Pikovsky, Physical Review E 102, 052315 (2020).
  * Clusella and Montbrió (2022) P. Clusella and E. Montbrió, arXiv preprint arXiv:2208.05515 (2022).
  * de Oliveira  _et al._ (2015) M. de Oliveira, M. da Luz, and C. Fiore, Physical Review E 92, 062126 (2015).
  * Luccioli  _et al._ (2012) S. Luccioli, S. Olmi, A. Politi, and A. Torcini, Physical review letters 109, 138103 (2012).
  * Restrepo  _et al._ (2005) J. G. Restrepo, E. Ott, and B. R. Hunt, Physical Review E 71, 036151 (2005).
  * Baer  _et al._ (1989) S. M. Baer, T. Erneux, and J. Rinzel, SIAM Journal on Applied mathematics 49, 55 (1989).
  * Ichinomiya (2004) T. Ichinomiya, Physical Review E 70, 026116 (2004).
  * Ko and Ermentrout (2008) T.-W. Ko and G. B. Ermentrout, Physical Review E 78, 016203 (2008).
  * Khasminskii (2011) R. Khasminskii,  _Stochastic stability of differential equations_ , Vol. 66 (Springer Science & Business Media, 2011).