# Stochastic synchronization induced by noise

Yunxiang Song<EMAIL_ADDRESS>Thomas A. Witten<EMAIL_ADDRESS>Department of Physics and James Franck Institute, University of Chicago,
Chicago, Illinois 60637, United States.

###### Abstract

Random perturbations applied in tandem to an ensemble of oscillating objects
can synchronize their motion. We study multiple copies of an arbitrary
dynamical system in a stable limit cycle, described via a standard phase
reduction picture. The copies differ only in their arbitrary phases $\phi$.
Weak, randomly-timed external impulses applied to all the copies can
synchronize these phases over time. Beyond a threshold strength there is no
such convergence to a common phase. Instead, the synchronization becomes
erratic: successive impulses produce stochastic fluctuations in the phase
distribution $q(\phi)$, ranging from near-perfect to near-random
synchronization. Here we show that the sampled entropies of these phase
distributions themselves form a steady-state ensemble, whose average can be
made arbitrarily negative by tuning the impulse strength. A random-walk
description of the entropy’s evolution accounts for the observed exponential
distribution of entropies and for the stochastic synchronization phenomenon.

##  I Introduction

Many kinds of signal transmission from cell phone and GPS signaling to
magnetic resonance imagingCallaghan (1991) depend on knowledge of the phase of
a remote oscillator. Because of this, it is valuable to find minimal ways to
control such phases. Remarkably, one can gain useful control by perturbing two
oscillators with random external forcing. In particular, an ensemble of
identical oscillators can be made to synchronize to a common phase by exposing
them all to the same random forcingPikovskii (1984); Pikovsky (1992); Pikovsky
and Kurths (1997); Teramae and Tanaka (2004); Goldobin and Pikovsky (2005);
Nakao  _et al._ (2007). Behavior consistent with this mechanism has been
reported in cortical neuronsMainen and Sejnowski (1995) and electronic
circuitsNagai and Nakao (2009). Here we examine a useful generalization of
this phenomenon in which the resulting synchronization is strong but
stochastic; the ensemble fluctuates between strong and weak synchronization.

We numerically analyzed this stochastic synchronization by repeatedly
disturbing the state of the oscillators using a specific common impulsive
force at random times, as in Ref Nakao  _et al._ (2005). The oscillators are
at different positions $\phi$ along their orbits at the moment of the impulse
and thus undergo different responses. Nevertheless, after a transient, each
oscillator returns to the limit-cycle orbit at some new position $\psi$
defined below. This $\psi$ depends on the orbit position $\phi$ at which the
impulse occurred; that is, $\psi$ is a function of $\phi$. For a given
oscillator and a given impulsive force, this “phase map” $\psi(\phi)$ is a
fixed function 111This phase map is denoted $F(\theta)$ in Ref. Nakao  _et
al._ (2005). This function suffices to determine the outcome after many
impulses. Many synchronization phenomena are amenable to this phase map or
“phase-reduction” analysis. One experimentally important type of forcing is
randomly-timed impulses with a common phase map, discussed in Sec II. An
initially random set of orbit positions $\phi$ thus repeatedly transforms into
a new set. Synchronization then amounts to progressive bunching of these
positions as the iteration proceeds. By contrast, the stochastic
synchronization treated here is more ambiguous.

Figure 1:  Polar plots of local density of phase angles during iteration using
a cubic phase map on the circle as described in Sec. III. The phase map used
has the form of Eq. (3) with $A=9.392$ and $\Phi=0.445$, for which the
Lyapunov exponent $\Lambda\simeq+0.141$ (cf. Sec. II). The initial ensemble
$\\{\phi^{1},...,\phi^{300}\\}_{0}$ consists of 300 positions sampled randomly
in an interval of width $0.00025$ at the upper left side of the circle. For
each sample $\alpha$, we denote the separation between $\phi^{\alpha}$ and its
successor around the circle as $\Delta\phi^{\alpha}$. We plot the $\log$ of
the sampled density i.e. $\log[1/(\Delta\phi^{\alpha})~{}]$ vs. the midpoint
position of the pair on the circle as a light-colored (red) dot. The
$\Delta\phi$’s are measured relative to that of uniformly-spaced $\phi$’s.
Circle radius corresponds to this uniform spacing. (a): Iterations 1, 21, 41,
…, 381 are shown by rows left to right. Radial spikes show regions of closely
spaced angles, contributing dominantly to the negative entropy $-H$ (Eq. (4)).
For comparison, dark colored (blue) dots show analogous $\phi$’s sampled from
a uniform distribution whose width is chosen to give the same entropy as the
phase-mapped $\phi$’s. For definiteness, these distributions are centered at
the upper left of the circle, like the initial distribution. (b): Detail of
iterations 201–205, showing progressive evolution over individual iterations.

Figure 1 illustrates the distinctive behavior of stochastic synchronization in
an ensemble of 300 oscillators subjected to a sequence of 380 randomly-timed
impulses, as described above. Each numbered circle depicts a set of $\phi$’s
at some iteration. Specifically, the light colored (red) dots around each
circle depict a sampling of the probability densities $q(\phi)$ inferred from
the spacings of adjacent $\phi$’s. Sharp spikes indicate phase angles having
high probability density. Thus synchronization of the ensemble would appear as
a single tall spike at each iteration, showing virtually the same phase angle
for all oscillators. However, in Fig. 1 there is no evident convergence to an
ordered, synchronized state. Instead, the spikes continually sharpen and
dissipate. Our aim in this paper is to quantify this ongoing stochastic
behavior of the ensemble of phases.

This behavior is distinct from the variety of synchronization phenomena
reported in the literature. It is quite distinct from the familiar
synchronization caused by mutual interaction of the oscillators, such as
mechanical clocks, fireflies, or applause Pikovsky  _et al._ (2003); Winfree
(1967); Mirollo and Strogatz (1990). Synchronization without such interaction
has been studied extensively by phase-reduction methodsTeramae and Tanaka
(2004); Nakao  _et al._ (2007). Noisy perturbations such as weak Gaussian
white noise or weak random telegraph noise was shown to lead to complete
synchronization. The breakdown of this synchronization for stronger
perturbations has been characterizedNagai  _et al._ (2005); Goldobin and
Pikovsky (2005); Pikovsky  _et al._ (2003). Our study below uses our methods
previously developed to study a colloidal realization of this transitionMoths
and Witten (2013); Eaton  _et al._ (2016); Witten and Diamant (2020).

In Sec. II to follow we provide a physical context for the random-timed
impulses and their statistical description. There we formally define our main
analytical tools: the Shannon entropy $H$ and average Lyapunov exponent
$\Lambda$. In Sec. III we describe the simulations that produced Fig. 1 in
terms of $H$ and $\Lambda$, noting that the distribution of $H$ values is a
broad exponential. In Sec. IV we give an argument to account for this observed
behavior of $H$ and predict how this distribution depends on $\Lambda$.
Lastly, in Sec. V, we relate stochastic synchronization to similar known
phenomena and suggest implications for future work.

##  II random-time impulses and cyclic motion

###  II.1 Phase coordinates

Our study considers physical systems that spontaneously repeat a periodic
cycle, such as a dripping faucet, a cooling iron, a beating heart, or a firing
neuronMainen and Sejnowski (1995). Such a system executes a periodic orbit in
its dynamical coordinate space in time with some period $T$. In stable
oscillators such as these, a moderate transient disturbance may move its
position away from the periodic orbit, but it eventually returns to the orbit.
We define a phase variable $\phi$ by first identifying a convenient starting
point on the orbit designated $\phi=0$. The unperturbed system advances from
this point with time $t$ and returns to the $\phi=0$ point in one period $T$.
We then label the orbit point reached at time $t$ by the “isochronous” phase
coordinate $\phi(t)\equiv t/T.$ In a random assembly of such oscillators all
times $t$ are equally represented. Therefore the probability distribution of
$\phi$, denoted $q(\phi)$ is uniform in $\phi$; i.e. $q(\phi)=1$.

We may determine the probability $q(\phi)$ by tracking a hypothetical ensemble
of many identical copies of the oscillator labeled by $\alpha$ . If the
resulting system is ergodic, this same $q(\phi)$ describes a single oscillator
over long times. Synchronization over time corresponds to a long-time bunching
of the various $\phi^{\alpha}$ towards a common value, which generally varies
in time. The corresponding $q(\phi)$ is sharply peaked. The unperturbed motion
simply advances all $\phi$’s around the phase circle at a uniform rate,
leaving their relative positions unchanged.

###  II.2 Lyapunov exponent and entropy

As noted above, the effect of a momentary disturbance on such an ensemble can
be summarized by a “phase map” $\psi(\phi)$ giving the relative phase long
after the disturbance, depending on the phase $\phi$ where the impulse
occurred, cf. Nakao  _et al._ (2005). Specifically. $\psi$ is the final phase
relative to that of a copy whose $\phi(t)=0$. It is convenient to designate
the point $\phi=0$ to be a fixed point where $\psi(\phi)=\phi$. Then a phase
map has the general appearance of Fig. 4, passing through the origin and
through the (equivalent) point (1, 1). Weak disturbances necessarily have
phase maps lying close to the $\psi=\phi$ line. No $\phi$ is changed very much
by the impulse, and any bunching effect is weak. The bunching effect on a
small interval of points near a given $\phi$ depends simply on the absolute
derivative $|\psi^{\prime}(\phi)|$ there. The corresponding interval of $\psi$
is widened by a factor $|\psi^{\prime}(\phi)|$. It is useful to define
$\lambda(\phi)$ to be the log of this widening factor. The average of this
$\lambda$ over all $\phi$ is denoted $\Lambda$ and is called the average
Lyapunov exponent
*[][.Our$Λ$isdenoted$λ$inthisreference.]nakao2005synchronyEaton  _et al._
(2016)

$\Lambda[\psi(\phi)]\equiv\int d\phi~{}(\lambda(\phi))=\int
d\phi~{}\left[\log\left|{d\psi\over d\phi}\right|\right]$ (1)

Weak impulses in general have $|\psi^{\prime}|$ near 1 and thus $\Lambda$
values near zero. The synchronization properties of a phase map depend
strongly on its $\Lambda$ as we explain below.

Further iterations of the same type result in the same mapping being applied
repeatedly to the ensemble. The result of such iteration is well documented in
textbooks on dynamical systems. The effects depend on delicate details such as
the precise amount of time elapsed between impulses. However, the effects on
disorder are more robust when the amount of time between impulses is random,
so that a given oscillator advances an arbitrary fraction of a cycle at the
next impulse. We term this process “random-shift iteration”. One may express
this operation as an alternating iteration of two maps. The first is the
$\psi(\phi)$ described above— identical for each iteration. The second is a
simple shift map $\tilde{\psi}(\phi)\equiv\phi+\beta_{i}$, the same for all
$N$ oscillators, but random for each iteration $i$. Thus the random shifts
$\beta_{i}$ produce a form of noise common to the oscillators. This random
shift limit corresponds to Poisson-distributed time between pulses in the
limit when their average rate goes to 0, as discussed in Ref. Nakao  _et al._
(2005).

To quantify the growth of order or disorder among the phase angles, we adopt
the entropy measure that gives the thermodynamic work needed to create the
order Reif (2009). The entropy $H$ of a probability measure $q(\phi)$ is
defined Shannon (1948) as

$H[q(\phi)]\equiv-\int d\phi~{}q(\phi)~{}\log q(\phi).$ (2)

This $H$ has a maximum value of 0 for uniform $q$, i.e. $q(\phi)=1$, with
complete disorder. Conversely, strongly bunched distributions have large,
negative $H$. One may estimate the entropy of an ensemble using a random
sampling of $\phi$ values from that ensemble, as detailed below.

###  II.3 Regimes of synchronization

Ref Eaton  _et al._ (2016) showed that for random-time iteration of a phase
map, the average change of entropy $\left\langle\Delta H\right\rangle$ in an
iteration is bounded above by the Lyapunov exponent $\Lambda$ of the map. For
strongly ordered ensembles with large negative $H$, these authors found that
$\left\langle\Delta H\right\rangle$ in an iteration approaches its upper bound
$\Lambda$. The study of Ref Eaton  _et al._ (2016) aimed to test these
predictions numerically using realistic phase maps calculated for a specific
soft-matter oscillator and a specific type of forcing, ranging from weak to
strong. As the amplitude of the forcing increased from zero, the $\Lambda$
decreased from zero, thus driving $H$ to arbitrarily negative values. Further
increase in the forcing gave a minimum $\Lambda$, where $H$ decreased most
rapidly. Still further increase in the forcing led to increasing $\Lambda$ and
slower decrease in $H$. Ultimately $\Lambda$ increased to positive values. For
small positive $\Lambda$ initial states with small $H$ initially showed the
predicted positive $\left\langle\Delta H\right\rangle$. But after many
iterations the ensemble reached a state of constant average $H$, denoted
$\left\langle H\right\rangle_{\infty}$. The authors noticed that this final
$\left\langle H\right\rangle$ could be made much smaller than the maximum
possible $H$ by making $\Lambda$ positive and close to zero. This is the
regime pictured in Fig. 1 that is the focus of the present work.

##  III Simulations

###  III.1 Methods

Our simulations aim to explore the transition between small negative
$\Lambda$, where $\left\langle H\right\rangle$ decreases indefinitely, and
small positive $\Lambda$, where $\left\langle H\right\rangle$ is nonzero but
seemingly arbitrarily small. We sought to verify that $\left\langle
H\right\rangle_{\infty}$ could be arbitrarily small yet nonzero. We sought to
determine the functional dependence of $\left\langle H\right\rangle$ on
$\Lambda$. And we sought to explore the generality of this behavior for
various map functions.

Since we aimed to explore general behavior, we used generic smooth phase map
functions rather than the experimentally motivated maps of Ref Eaton  _et al._
(2016). Accordingly, we studied a qualitatively similar cubic function of the
form

$\psi(\phi)\equiv\phi+A\phi(1-\phi)(\Phi-\phi).$ (3)

Here the amplitude parameter $A$ regulates the distance between $\phi$ and
$\psi$, and $\phi=\Phi$ is a fixed point where $\psi(\phi)=\phi$. By
increasing the amplitude $A$ one could increase the value of $\Lambda$. All
$\psi$’s differing by an integer denote the same point on the circle.

Given this map function, we formally describe the step-by-step procedure for
carrying out our numerical simulations. This procedure allows us to judge the
dependence on amount of iteration, and its variation with the sequence of
random inputs $\beta_{i}$.

  1. 1.

Ensemble initialization: An ensemble of $N$ oscillator phases
$\\{\phi^{\alpha}\\}_{0}$ where $\alpha=1,2,\dots N$ is drawn from a uniform
probability density $q_{0}(\phi)$ supported on the interval
$[\phi_{0},\phi_{0}+w)$ on the unit circle of width $w$. The left endpoint
$\phi_{0}$ of this interval is randomly selected.

  2. 2.

Ensemble iteration: Each phase $\phi^{\alpha}$ in the ensemble is mapped to an
image phase under the phase map $\psi(\phi)$ (Eq. (3)). A random shift
$\beta_{0}$ is applied to all phase in the ensemble to obtain
$\phi^{\alpha}_{1}\equiv\psi(\phi_{0}^{\alpha})+\beta_{0}$, with a
corresponding probability density $q_{1}(\phi)$.

  3. 3.

Entropy estimation: The entropy $H$ is estimated from the discrete collection
of phases $\\{\phi^{\alpha}\\}$ Victor (2002). When the number of sampled
positions $N$ is large, this estimate can be written

$H\left[\\{\phi^{\alpha}\\}\right]\simeq\frac{1}{N}\sum_{\alpha=1}^{N}\log\delta^{\alpha}+\log\left[2N-2\right]+\gamma~{},$
(4)

where $\delta^{\alpha}$ is the nearest-neighbor distance for phases
$\phi^{\alpha}$ and $\gamma\simeq 0.577$ is the Euler-Mascheroni constant.

  4. 4.

Repeated iteration: steps 2 and 3 are repeated for $n$ iterations, with a
random sequence of shifts $\beta_{i}$, to generate a single sequence of
$\\{\phi^{\alpha}\\}_{i}$ and corresponding entropies $H_{i}$, for i=1,2,
…,$n$. We typically choose $n$=200.

  5. 5.

Repeated sequence: step 4 is repeated for $k$ randomly generated sets of $n$
iterations, to obtain $k$ independent trajectories. We typically choose
$k$=100.

  6. 6.

Averaging over sequences: The entropy at each iteration averaged over all $k$
sets of iterations gives a trajectory for the average entropy $\left\langle
H\right\rangle_{i}$ as a function of iteration $i$.

###  III.2 Results

In agreement with Eaton  _et al._ (2016) we observe three types of
characteristic behavior for iterated phase maps $\psi(\phi)$ with different
$\Lambda$. We plot examples of these distinct behaviors in Fig. 2.

Figure 2:  Entropy $H$ vs iteration number $i$ for randomly-shifted iterated
circle maps as described in Sec. III.1. Panels (a)–(c) were made from maps of
the form of Eq. (3) with three different values of Lyapunov exponent $\Lambda$
as indicated. Each color tracks the entropy change over a single trajectory.
Black trajectory shows the entropy at each iteration averaged over all colored
trajectories. Straight colored line at left indicates slope $\Lambda$.
Concentration of points for $H\simeq-33$ is due to limited machine precision,
as described in the text.

The first class of behavior is for negative $\Lambda$, as shown in Fig. 2(a).
Here the average entropy decreases to arbitrarily small values over time, so
that the system becomes synchronized. The $\left\langle H\right\rangle$
decreases at a rate consistent with $\Lambda$, as explained e.g. in Ref. Eaton
_et al._ (2016). Here, we note the entropy floor in Figs. 2a,b just below
$H=-30$ instead of an indefinite decay towards $-\infty$. Such a floor is an
artifact of our simulation implementation. When two neighboring phase angles
are indistinguishable up to machine precision, we manually separate them by
1$\times 10^{-15}$ to maintain finite values in the the sum of Eq. (4). The
floor is reached when all phase angles in the ensemble are indistinguishable.
These numerical limitations produced significant distortions in our measured
$\left\langle H\right\rangle_{\infty}$ for $\Lambda\mathrel{\raise
1.29167pt\hbox{$<$\kern-7.5pt\lower 4.30554pt\hbox{$\sim$}}}0.14$ (cf. Fig.3).

The second class of behavior is for large positive $\Lambda$, as shown in Fig.
2(c). In this case, the average entropy approaches the maximally disordered
state with $\left\langle H\right\rangle\simeq 0$. Such behavior is also well
understood Eaton  _et al._ (2016); Witten and Diamant (2020).

The third class of behavior exists in a narrow regime between the previous two
classes, where $\Lambda$ is small, but still positive, so that there is an
initial upwards slope in $\left\langle H\right\rangle$, as shown in Fig. 2(b).
The average change in entropy initially follows the upper bound $\Lambda$ but
eventually reaches a steady state. The small asymptotic value $\simeq-5$
indicates states that are on average much more ordered than a completely
random set of phases.

We made further statistical measurements of the entropy in this steady-state
regime of small average entropy. First, we verified that $\left\langle
H\right\rangle_{\infty}$ indeed extrapolates to $-\infty$ as
$\Lambda\rightarrow 0$ in the undistorted range $\Lambda>0.14$ (Fig. 3). The
limiting behavior was consistent with $\left\langle
H\right\rangle_{\infty}\sim 1/\Lambda$. Changing the map while keeping
$\Lambda$ fixed had little effect on $\left\langle H\right\rangle_{\infty}$,
as shown in Fig. 4. Next, we measured the distribution of entropy values
$P(H)$, shown in Fig. 5. The measurements strongly indicate an exponential
falloff with $|H|$. Finally, we measured the statistics of the incremental
change $h_{i}\equiv H_{i+1}-H_{i}$, as shown in Fig. 6. For the most negative
range of $H$, the mean $\left\langle h\right\rangle$ of this distribution
approached $\Lambda$, as anticipated Eaton  _et al._ (2016), while the
variance was of order unity. The $h$ values fell off sharply on both sides of
the mean, with no apparent long tails.

Figure 3:  Steady-state entropy vs. average Lyapunov exponent. Dark color
(blue) dots plot simulated data for a selection of iterated maps with near
continuous range of $\Lambda$. Light color (red) curve is a fit to the form
$\left\langle H\right\rangle_{\infty}=a/\Lambda$ anticipated in Sec. IV with
$a=0.77$. Data for small $\Lambda$ was subject to distortion in the positive
direction owing to our treatment of small separations of the
$\\{\phi^{\alpha}\\}$ cf. Sec. III.2.  Figure 4:  Phase maps $\psi(\phi)$ of
the form of Eq. (3), used to study the effect of varying the phase map with
fixed $\Lambda\simeq 0.141$. Fixing $\Lambda$ allows one free parameter within
the class of cubic maps. The $\left\langle H\right\rangle_{\infty}$ values
obtained from simulating 30 maps spanning this range was the same within the
statistical uncertainty, viz. $\left\langle H\right\rangle_{\infty}=-5.5\pm
0.2$ . Figure 5:  Cumulative distribution of $H$ as obtained in Sec. III under
conditions of Fig.2b, using a map with Lyapunov coefficient $\Lambda=0.14$,
shown as a semilogarithmic plot. The $P(H)$ defined in the text is the
derivative of the cumulative distribution plotted. Dark line shows for each
entropy $H$ the fraction of samples with entropy smaller than $H$. Mid section
of this curve is approximately a straight line indicating exponential falloff.
Light line is the exponential distribution with a scale height $H_{0}$
matching this section: $H_{0}=3.85$ . Figure 6:  Histograms of the
distribution of incremental changes of $H$, $p(h)$ for the map with
$\Lambda=0.14$. Light-colored histogram includes the full sample. Middle color
and foreground color show the half of the sample with smallest $H$ and the 20
percent of the sample with the smallest $H$. The variances of all the samples
were approximately unity. The mean for the whole sample was approximately 0,
as expected for a steady-state process. For the foreground sample, the mean
$\left\langle h\right\rangle\simeq 0.16$, is equal to $\Lambda$ within
sampling uncertainty. The two sharp dropoffs near $h=0.4$ and $h=1.8$
correspond to extremal slopes in $\psi(\phi)$.

##  IV Explanation of $H$ Distribution in steady state

Here we show that the regime of stochastic fluctuations of the phases with
indefinitely small entropy is to be expected for general random-shift phase
maps with small positive $\Lambda$. We present a scheme that accounts for this
behavior. It accounts for the observed broad distribution of $H$ values in the
steady state. It also predicts that the average $\left\langle
H\right\rangle_{\infty}$ varies as $1/\Lambda$ for $\Lambda\rightarrow 0$.

Both the asymptotic scaling of $\left\langle H\right\rangle$ for small
Lyapunov exponent $\Lambda$ and the functional form of the distribution $P(H)$
can be understood in stochastic terms: the fluctuations of $H_{i}$ may be
viewed as a random walk process. In this regime, the typical $H$ values and
their average are indefinitely smaller than the values observed for large
$\Lambda$. Our description exploits the following features of this
small-$\Lambda$ regime:

  1. 1.

The incremental changes in $H$ in one iteration $i$, $h_{i}\equiv
H_{i+1}-H_{i}$, vary significantly, but the range of this variation is limited
and is little dependent on $H_{i}$ when $H_{i}\ll 0$.

  2. 2.

The mean value of the increment $h$ approaches $\Lambda$ as
$\Lambda\rightarrow 0$, as shown in Ref Eaton  _et al._ (2016).

  3. 3.

$H$ can never exceed a maximum, namely its value for a uniform probability
distribution of phase angles. With our conventions this maximum is 0.

These conditions suggest a simple diffusion-drift mechanism Einstein (1905);
Weiss and Rubin (1982) for the strong, steady-state entropy fluctuations that
we observe. We picture the changes of $H$ as a slightly biased random walk,
whose steps $h_{i}$ are drawn independently from a distribution $p(h)$. This
$p(h)$ is taken to be independent of $H$, as motivated by feature 1) above.
The mean of this distribution $\left\langle h\right\rangle$ approaches
$\Lambda$, as dictated by feature 2. Finally, any step that leads to $H>0$ is
withdrawn, as dictated by feature 3. In the absence of the drift, the entropy
undergoes unbiassed excursions and ultimately reaches indefinitely negative
values. However, with small positive $\left\langle h\right\rangle$, these
excursions are opposed by the drift. Then $H$ can no longer decrease
indefinitely, and a steady state is reached in which the drift is balanced by
the random steps. This confinement is weaker and $\left\langle H\right\rangle$
is more negative as $\left\langle h\right\rangle\rightarrow 0$. Thus in the
small-positive-$\Lambda$ regime of interest, $\Lambda^{2}\ll\left\langle
h^{2}\right\rangle\ll\left\langle H\right\rangle_{\infty}^{2}$.

The behavior of such diffusion-drift processes is well known Weiss and Rubin
(1982) . In particular, in steady state, the probability distribution $P(H)$
obeys

${d^{2}P\over d~{}H^{2}}={2\left\langle h\right\rangle\over\left\langle
h^{2}\right\rangle}~{}{dP\over dH}~{}~{},$ (5)

so that

$P(H)=C~{}\exp{(H/H_{0})}~{}~{},$ (6)

where $C$ is a normalizing constant and the scale height $H_{0}$ is given by

$H_{0}=\left\langle h^{2}\right\rangle/(2\left\langle h\right\rangle).$ (7)

This predicted behavior is consistent with the simulations described above.
Fig. 5 shows the observed distribution of $H$ values for $\Lambda=0.14$ from
Sec. III, incorporating 300 runs of 200 time steps apiece. Time steps before
step 100 were discarded to avoid effects of the initial transient; the
remaining sample contained $10^{4}$ measurements. The plotted distribution is
the fraction of samples with entropy more negative than $H$. The picture shows
a one-decade range consistent with an exponential distribution as anticipated
in Eq. 6, implying $H_{0}\simeq 3.85$ .

For a more explicit test of the explanation, we measured the distribution of
step changes $h$ observed in the simulation. In general this distribution
depends on $H$. Since good prediction is only expected for
$H\ll-\sqrt{\left\langle h^{2}\right\rangle}$, we measured incremental $h$
values only for the smallest 20 percent of the $H$ values, as shown in Fig. 6.
Here we obtained $\left\langle h\right\rangle\simeq 0.16$. This value violates
the upper bound Eaton  _et al._ (2016) on the true average relative to the
Lyapunov coefficient $\Lambda$ of the map defined in Eq. (1), viz. ,
$\left\langle h\right\rangle<\Lambda=0.14$. We attribute this discrepancy to
statistical uncertainty in the simulation 222 According to the diffusion-drift
picture, we expect an uncertainty in $\left\langle h\right\rangle$ of order
$\left\langle h^{2}\right\rangle/\sqrt{N}$, where $N(=2000)$ is the number of
samples. Thus uncertainty is roughly 0.02, comparable to the discrepancy
between 0.16 and 0.14. . Thus, we regard the measured $\left\langle
h\right\rangle$ as a crude consistency check. The variance $\left\langle
h^{2}\right\rangle$ was approximately $1.0$. These figures give an expected
$H_{0}$ using Eq. 7: $H_{0}=3.1$ using the measured $\left\langle
h\right\rangle$, or $H_{0}=3.5$ using $\left\langle h\right\rangle=\Lambda$.
This is roughly consistent with the observed $H_{0}=3.85$ shown in Fig 5, thus
lending support to the diffusion-drift picture.

The diffusion-drift picture also accounts for the dependence of $\left\langle
H\right\rangle$ on $\Lambda$ reported in Sec. III. For asymptotically small
$\Lambda$, we may replace $\left\langle h\right\rangle$ by $\Lambda$. The
diffusion-drift picture then gives an exponential distribution $P(H)$ with
scale height $H_{0}\rightarrow\left\langle h^{2}\right\rangle/(2\Lambda)$ and
a mean value $\left\langle H\right\rangle\rightarrow-H_{0}$, which varies as
$\Lambda^{-1}$, in agreement with Fig. 3.

Our actual $\Lambda$ values did not reach this asymptotic regime. Even our
smallest $\Lambda$ simulation had a $P(H)$ that departed significantly from an
exponential as seen in Fig. 5. It had a mean $\left\langle
H\right\rangle=-4.8$ that was somewhat smaller than $-H_{0}~{}(=-3.85)$. The
predicted $\left\langle H\right\rangle\rightarrow\left\langle
h^{2}\right\rangle/(2\Lambda)$ is qualitatively consistent with the
measurements of Fig 3. We verified that $\left\langle h^{2}\right\rangle\simeq
1$ for large negative $H$ over the range of $\Lambda$ studied. Thus the
asymptotic formula predicts $\left\langle H\right\rangle$ of the form
$a/\Lambda$, in agreement with the fitted curve in Fig. 3. However, the
predicted value of the coefficient $a$ is 0.5—only 64 percent of the fitted
value of $a$.

Overall, the asymptotic diffusion-drift picture appears to account adequately
for the limited simulations reported above. In the next section we discuss the
generality of this picture.

##  V Discussion

In the foregoing we have demonstrated a novel mechanism for synchronizing an
ensemble of oscillators with random phases using a robust, physically
applicable forcing protocol. The procedure can produce controlled disorder of
nonzero magnitude as well as indefinitely small disorder. In this discussion
section we justify why our low-entropy states can be viewed as a form of
synchronization. We argue that our system achieves a statistically well-
behaved though unusual steady state, and cite a similar established example.
Finally we review the practical limitations of the method, and survey possible
generalizations.

###  V.1 Entropy and synchronization

Synchronization is conventionally quantified by the uncertainty spread in
between the phase angles of two oscillators. However a single narrow interval
of uncertainty is not necessary in order to have useful knowledge about a
phase. If the probability is divided into two narrow intervals, the phase
information does not degrade greatly simply by separating the two intervals by
a large spacing. Indeed, the number of trials needed to determine the phase
angle to a given tolerance is equal for the single interval and the split
interval regardless of the size of the split. The difference in entropy of two
distributions is precisely the relative number of trials needed on average for
the two. Thus for any given distribution $q(\phi)$ one may identify a uniform
distribution that is equivalent in number of trials needed to determine the
phase. This is the uniform distribution whose entropy is equal to that of
$q(\phi)$. For many purposes of information transmission, this entropy measure
is what determines the transmitting capacity Shannon (1948).

###  V.2 Statistical regularity

The erratic variability of the synchronization and of the entropy in our
system raises the question whether the process is even statistically well-
defined. Several tests gave us positive reassurance on this point. The
behavior of the simulated entropy fluctuations appears consistent with a
stationary distribution $P(H)$. On the one hand, the measured average entropy
showed a variance that decreases inversely with number of iterations or
trials, as expected for averaging from a fixed ensemble.

On the other hand the distribution of $\lambda(\phi)$
($\equiv\log|\psi^{\prime}(\phi)|$) values for generic functions $\psi(\phi)$
is also consistent with a well-defined $P(H)$ distribution. When $H_{i}$ is
small, this $\lambda(\phi)$ determines the distribution of $H_{i+1}$ values in
one iteration starting from a given distribution $q_{i}(\phi)$. The particular
value of $H_{i+1}$ obtained depends on the randomly chosen shift amount
$\beta_{i}$ (Sec. II.2). This dependence can be understood in the limit where
$q_{i}(\phi)$ has support confined to a few narrow intervals on the circle.
Then the $H_{i+1}(\beta_{i})$ for all $\beta_{i}$ can be readily inferred from
the properties of $\lambda(\phi)$, as shown in the Appendix. In the opposite
regime of uniform $q_{i}$, there is no $\beta$ dependence of $H_{i+1}$ and
thus no randomness in the distribution of $h$. Thus the assumption of well-
behaved randomness of $h$ is confirmed for the highest as well as the lowest
initial entropy $H_{i}$.

###  V.3 Parallels

Stochastic synchronization extends the notion of “noisy on-off intermittency”
or “bursting behavior” reported in the noise-induced synchronization
literature Teramae and Tanaka (2004); Nakao  _et al._ (2005, 2007). This
behavior was seen in oscillators subject to two types of random perturbation:
(a) common impulses as described above with negative Lyapunov exponent and (b)
uncorrelated Langevin noise. There the intermittency arose from the
competition between the disruptive influence of the noise and the
synchronizing effect caused by the local Lyapunov exponent $\lambda(\phi)$
over many iterations. Our work shows that this intermittent steady state
extends to positive $\Lambda$, where the disordering effect of $\Lambda$ is
opposed by the ordering effect of any phase map on a near-uniform distribution
$q(\phi)$.

More broadly, bursting intermittency is encountered in dynamical systems near
a transition to chaos. Properties of this intermittency such as scaling of
correlations and size of the chaotic region have been explained using
properties of the incipient chaotic state Heagy  _et al._ (1994);
Venkataramani  _et al._ (1996); Yu  _et al._ (1990). The present work doesn’t
require proximity to a chaotic state and uses only the phase map associated
with some perturbation. Nevertheless, the nature of our disorder appears
similar to that seen in the prior works. It seems likely that the degree of
disorder in these near-chaotic systems quantitatively resembles that shown in
the present system, that operates well within the boundaries of a stable limit
cycle.

Stochastic states analogous to ours are well known in disordered wave systems.
There, as in our system, configurational variables analogous to $\phi$ lack
well-defined probabilities that reveal the statistical behavior of the system.
One case of similar behavior is the transmission of a wave such as light
through a stack of $n$ different transmitting layers such as sheets of glass
of different thicknesses Stone and Joannopoulos (1981); Chabanov  _et al._
(2000). If the sheets are all much thicker than a wavelength, the complex
reflection and transmission amplitudes vary widely from sheet to sheet. As a
result, the overall transmission varies erratically with wavelength, with no
single characteristic distribution. The distribution of transmission
coefficients is ill-conditioned and only the logarithm of the transmission
coefficient has well-conditioned statistics. In our system the probability
distribution $q(\phi)$ is ill-conditioned and the entropy does not reach a
definite value for large systems. Instead the entropy converges to a steady-
state distribution with a well-determined average.

###  V.4 Limitations

In investigating this class of stochastic dynamics we encountered limitations
in exploring the asymptotic regime of strong but incomplete ordering.
Numerically, we were unable to reliably explore maps with positive Lyapunov
exponents $\Lambda$ much smaller than 0.14 and corresponding $\left\langle
H\right\rangle\mathrel{\raise 1.29167pt\hbox{$<$\kern-7.5pt\lower
4.30554pt\hbox{$\sim$}}}-5$. The difficulty is already apparent in Figs. 2a,b.
There, as noted in Sec. III.2 one sees several trajectories that descend to
the bottom of the figure and remain there. These trajectories have
$\phi^{\alpha}$ values separated by less than the machine precision of the
calculations. Such pairs of $\phi$’s cannot give valid contributions to Eq.
(4). These trajectories necessarily increase in number as the iterations
proceed, eventually compromising the measurement of entropy. Thus the number
of iterations attainable in practice is limited. This in turn limits the
$\Lambda$ values that can be explored, since small $\Lambda$ entails slow
relaxation to the steady state. These numerical limitations can be mitigated
by increasing the numerical precision of the simulation. In this way one could
improve our crude validation of the mechanism described in Sec. IV . This
improvement would be very desirable.

Similar limitations would be expected in experiments like the soft-matter
experiments treated in Ref. Eaton  _et al._ (2016). First, the final state of
any dynamical oscillator is subject to random noise as well as to the imposed
random-time phase maps. This noise requires the forcing to be sufficiently
strong to induce rapid synchronization. Any synchronization is thus limited by
the random noise Nakao  _et al._ (2007). A second limitation comes from the
inevitable variability of colloids or other oscillators in the ensemble. This
variability means each element $\alpha$ of the ensemble has a distinct phase
map $\psi_{\alpha}(\phi)$.

Such limitations mean that clear-cut experimental attainment of the asymptotic
predictions above is unlikely. Nevertheless, this limit could be useful for
devising well-behaved statistics to describe the stochastic synchronization
regime, as was the case with disordered quantum wires in Ref. Stone and
Joannopoulos (1981).

###  V.5 Generalizations

The stochastic dynamics treated here was chosen for its simplicity. It treats
only a very restricted type of disruption of the limit cycle. Still, the
anomalous ordering behavior demonstrated here should occur more generally.
Many of the known noise-induced synchronization phenomenaTeramae and Tanaka
(2004); Nakao  _et al._ (2005); Yoshida  _et al._ (2006); Nakao  _et al._
(2007); Hata  _et al._ (2010) have been demonstrated for broad classes of
noise, including continuous random driving. These phenomena are regulated by
an average Lyapunov exponent $\Lambda$ analogous to ours. Generally
synchronization was reported only when this $\Lambda$ was negative. Our
results suggest that a form of stochastic synchronization can occur even when
$\Lambda$ is positive.

For some incremental generalizations of our simple forcing, stochastic
synchronization like that shown here appears likely. For example, if the
impulsive forces are not restricted to be identical, but instead vary in
strength, the ordering should persist. Such variable forcing would lead to a
range of phase maps $\psi_{i}(\phi)$ for different iterations $i$. Still if
these maps all had $\Lambda$ values bounded above by a small, positive number,
one would expect bounds on the entropy similar to those shown here. Likewise,
we restricted the impulses to be widely separated in time so that the
oscillators all relax to their periodic cycles before the next impulse. This
requirement simplifies our analysis but does not appear essential to the
phenomenon of stochastic ordering Nakao  _et al._ (2005).

One general virtue of noise-induced synchronization is that one may use it to
induce synchronization without specific knowledge about the limit cycle being
synchronized, such as the phase map function or the cycle time. We note that
this virtue is preserved for the impulsive noise of the current study.

##  VI Conclusion

As noted in the Introduction, synchronization of a remote oscillator enables
transmission of information. Now, our demonstration of low-entropy oscillator
states is far from a demonstration of effective information transmission. Yet
this stochastic counterpart of synchronization shows potential as a
generalized means of transmission, potentially useable by technology or living
systems. Further, this model of intermittency may offer a tool to address open
questions in wave localization and strong turbulence.

## Acknowledgments

We thank Ivar Martin, Peter Littlewood, Kyle Kawagoe, and Alex Edelman for
insightful discussions. Jonah Eaton and Martin Falk suggested valuable
improvements in the manuscript. Y.S. acknowledges support from the University
of Chicago Jeff Metcalf Internship Program.

## Appendix: Limits on the range of $H_{i+1}$ values

Here we address the claim made in Sec. V.2 that the change in entropy under a
single iteration of our random-time phase mapping varies within a limited
range. We consider the regime where the initial distribution $q_{i}(\phi)$ is
narrowly defined, as is typical when the initial entropy $H_{i}$ is small. In
general the new probability distribution $q_{i+1}(\hat{\phi})$ can be
expressed Eaton  _et al._ (2016) as

$q_{i+1}(\hat{\phi})=\int
d\phi~{}q_{i}(\phi)~{}\delta\left(\hat{\phi}-\psi(\phi+\beta_{i})\right).$ (8)

Then $H_{i+1}$ is obtained from $q_{i+1}(\hat{\phi})$ by applying Eq. (2). We
note that when $\psi(\phi)$ is non-monotonic, several values of $\phi$ may
contribute to a given $\hat{\phi}$.

As announced in the main text, we first consider the range of $H_{i+1}$ values
when $q_{i}$ is confined to an arbitrarily small support. For such “well-
ordered” $q_{i}(\phi)$ there is generally only one $\phi$ contributing to each
$\hat{\phi}$ in Eq. (8). Then Eq. (8) becomes

$\displaystyle q_{i+1}(\hat{\phi})=q_{i}(\phi)$ $\displaystyle~{}[\int
d\phi~{}\delta(~{}\hat{\phi}-\psi(\phi+\beta_{i})~{})]$
$\displaystyle=q_{i}(\phi)$
$\displaystyle~{}|\psi^{\prime}(\phi+\beta_{i})|^{-1}.$ (9)

For such $q_{i+1}$, $H_{i+1}(\beta_{i})$ becomes a convolution of
$\lambda(\phi+\beta)$ with $q_{i}(\phi)$. To show thisEaton  _et al._ (2016),
we use the definition from Eq. (2) with Eq. (Appendix: Limits on the range of
$H_{i+1}$ values) and the fact that
$d\hat{\phi}~{}q_{i+1}(\hat{\phi})=d\phi~{}q_{i}(\phi)$ to write $H_{i+1}$ as

$H_{i+1}(\beta_{i})=-\int
d\phi~{}q_{i}(\phi)~{}\log\left(~{}|\psi^{\prime}(\phi+\beta_{i})|^{-1}~{}q_{i}(\phi)~{}\right).$
(10)

After decomposing the $\log$ and noting that
$\log(|\psi^{\prime}(\phi+\beta_{i})|)\equiv\lambda(\phi+\beta_{i})$, this
gives

$H_{i+1}(\beta_{i})=\int d\phi~{}q_{i}(\phi)~{}\lambda(\phi+\beta_{i})-\int
d\phi~{}q_{i}(\phi)~{}\log q_{i}(\phi).$ (11)

Here the second term (simply $H_{i}$) is independent of $\beta_{i}$ while the
first term is the convolution claimed.

We now argue that the range of $H_{i+1}(\beta_{i})$ values is strongly
limited. This range is greatest when $q_{i}(\phi)$ is simply a delta function.
Then, ignoring the constant second term in Eq. (11), the distribution of
$H_{i+1}$ values is simply the distribution of $\lambda$ values. This
distribution may be readily evaluated for any given phase map $\psi(\phi)$.
For a smooth map function like those encountered in Eaton  _et al._ (2016) and
those considered here, there is a maximum $\lambda$ at the point of maximum
absolute slope. There is no minimum value since there are extrema where $\phi$
has arbitrarily small $|\psi^{\prime}|$. We consider generic $\psi(\phi)$ with
no points where $\psi^{\prime}$ and $\psi^{\prime\prime}$ simultaneously
vanish, so that all extrema of $\psi(\phi)$ are quadratic. Near any extremum
$\phi^{*}$, $|\psi^{\prime}|\sim|\phi-\phi^{*}|$, so that
$\lambda(\phi)\sim\log|\phi-\phi^{*}|$. Then the $p(\lambda)$ falls off
exponentially for large negative $\lambda$. Thus for well-ordered
$q_{i}(\phi)$ the range of $H_{i+1}$ is well-confined, the central limit
theoremWeiss and Rubin (1982) applies, our diffusion hypothesis of Sec. IV is
well justified and a well-defined $P(H)$ is reached.

The range of $H_{i+1}(\beta_{i})$ is also narrow in the opposite limit of
uniform $q_{i}(\phi)$. Here the shift $\beta_{i}$ has no effect on
$q_{i}(\phi)$ and thus $H_{i+1}(\beta_{i})$ is independent of $\beta_{i}$; The
spread of $H_{i+1}$ values for a given $H_{i}$ goes to zero. While these
arguments cover only the limiting cases and fall short of a proof, they lend
plausibility to our numerical finding that the distribution of entropy
increments $h_{i}$ is statistically well-behaved.

## References

  * Callaghan (1991) Paul T. Callaghan,  _Principles of nuclear magnetic resonance microscopy_ (Clarendon Press; Oxford University Press, New York,, 1991).
  * Pikovskii (1984) A. S. Pikovskii, “Synchronization and Stochastization of Nonlinear Oscillations by External Noise,” in  _Nonlinear and Turbulent Processes in Physics_ (1984) p. 1601.
  * Pikovsky (1992) Arkady S. Pikovsky, “Statistics of trajectory separation in noisy dynamical systems,” Physics Letters A 165, 33–36 (1992).
  * Pikovsky and Kurths (1997) Arkady S. Pikovsky and Jürgen Kurths, “Coherence resonance in a noise-driven excitable system,” Phys. Rev. Lett. 78, 775–778 (1997).
  * Teramae and Tanaka (2004) Jun-nosuke Teramae and Dan Tanaka, “Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators,” Phys. Rev. Lett. 93, 204103 (2004).
  * Goldobin and Pikovsky (2005) Denis S. Goldobin and Arkady Pikovsky, “Synchronization and desynchronization of self-sustained oscillators by common noise,” Phys. Rev. E 71, 045201 (2005).
  * Nakao  _et al._ (2007) Hiroya Nakao, Kensuke Arai, and Yoji Kawamura, “Noise-induced synchronization and clustering in ensembles of uncoupled limit-cycle oscillators,” Physical review letters 98, 184101 (2007).
  * Mainen and Sejnowski (1995) Zachary F Mainen and Terrence J Sejnowski, “Reliability of spike timing in neocortical neurons,” Science 268, 1503–1506 (1995).
  * Nagai and Nakao (2009) Ken Nagai and Hiroya Nakao, “Experimental synchronization of circuit oscillations induced by common telegraph noise,” Physical Review E 79, 036205– (2009).
  * Nakao  _et al._ (2005) Hiroya Nakao, Ken-suke Arai, Ken Nagai, Yasuhiro Tsubo, and Yoshiki Kuramoto, “Synchrony of limit-cycle oscillators induced by random external impulses,” Physical Review E 72 (2005), 10.1103/physreve.72.026220.
  * Pikovsky  _et al._ (2003) Arkady Pikovsky, Jurgen Kurths, Michael Rosenblum, and Jürgen Kurths,  _Synchronization: a universal concept in nonlinear sciences_ , Vol. 12 (Cambridge university press, 2003).
  * Winfree (1967) Arthur T Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” Journal of theoretical biology 16, 15–42 (1967).
  * Mirollo and Strogatz (1990) Renato E Mirollo and Steven H Strogatz, “Synchronization of pulse-coupled biological oscillators,” SIAM Journal on Applied Mathematics 50, 1645–1662 (1990).
  * Nagai  _et al._ (2005) Ken Nagai, Hiroya Nakao, and Yasuhiro Tsubo, “Synchrony of neural oscillators induced by random telegraphic currents,” Physical Review E 71 (2005), 10.1103/physreve.71.036217.
  * Moths and Witten (2013) Brian Moths and TA Witten, “Full alignment of colloidal objects by programed forcing,” Physical review letters 110, 028301 (2013).
  * Eaton  _et al._ (2016) Jonah A Eaton, Brian Moths, and Thomas A Witten, “Criterion for noise-induced synchronization: application to colloidal alignment,” Physical Review E 94, 032207 (2016).
  * Witten and Diamant (2020) Thomas A Witten and Haim Diamant, “A review of shaped colloidal particles in fluids: anisotropy and chirality,” Reports on Progress in Physics 83, 116601 (2020).
  * Reif (2009) Frederick Reif,  _Fundamentals of statistical and thermal physics_ (Waveland Press, 2009).
  * Shannon (1948) C E Shannon,  _The mathematical theory of communication._ (University of Illinois Press, 1948).
  * Victor (2002) Jonathan D Victor, “Binless strategies for estimation of information from neural data,” Phys. Rev. E 66, 051903 (2002).
  * Einstein (1905) Albert Einstein, “On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat,” Annalen der physik 17, 549–560 (1905).
  * Weiss and Rubin (1982) George H. Weiss and Robert J. Rubin, “Random walks: Theory and selected applications,” in  _Advances in Chemical Physics_ (John Wiley & Sons, Ltd, 1982) pp. 363–505.
  * Heagy  _et al._ (1994) J. F. Heagy, N. Platt, and S. M. Hammel, “Characterization of on-off intermittency,” Phys. Rev. E 49, 1140–1150 (1994).
  * Venkataramani  _et al._ (1996) Shankar C. Venkataramani, Thomas M. Antonsen, Edward Ott, and John C. Sommerer, “On-off intermittency: Power spectrum and fractal properties of time series,” Physica D: Nonlinear Phenomena 96, 66–99 (1996), measures of Spatio-Temporal Dynamics.
  * Yu  _et al._ (1990) Lei Yu, Edward Ott, and Qi Chen, “Transition to chaos for random dynamical systems,” Phys. Rev. Lett. 65, 2935–2938 (1990).
  * Stone and Joannopoulos (1981) A. Douglas Stone and J. D. Joannopoulos, “Probability distribution and new scaling law for the resistance of a one-dimensional anderson model,” Phys. Rev. B 24, 3592–3595 (1981).
  * Chabanov  _et al._ (2000) A. A. Chabanov, M. Stoytchev, and A. Z. Genack, “Statistical signatures of photon localization,” Nature 404, 850–853 (2000).
  * Yoshida  _et al._ (2006) Katsutoshi Yoshida, Keijin Sato, and Atsushi Sugamata, “Noise-induced synchronization of uncoupled nonlinear systems,” Journal of sound and vibration 290, 34–47 (2006).
  * Hata  _et al._ (2010) Shigefumi Hata, Takeaki Shimokawa, Kensuke Arai, and Hiroya Nakao, “Synchronization of uncoupled oscillators by common gamma impulses: from phase locking to noise-induced synchronization,” Physical Review E 82, 036206 (2010).