# Extreme multistability in symmetrically coupled clocks

Zhen Su  Potsdam Institute for Climate Impact Research, 14473 Potsdam, Germany
Department of Computer Science, Humboldt-Universität zu Berlin, 10099 Berlin,
Germany  Jürgen Kurths  Potsdam Institute for Climate Impact Research, 14473
Potsdam, Germany  Department of Physics, Humboldt-Universität zu Berlin, 10099
Berlin, Germany  Yaru Liu<EMAIL_ADDRESS>Potsdam Institute for Climate
Impact Research, 14473 Potsdam, Germany  Department of Mathematics, Jinan
University, 510632 Guangzhou, China  Serhiy Yanchuk  Potsdam Institute for
Climate Impact Research, 14473 Potsdam, Germany  Institute of Mathematics,
Humboldt-Universität zu Berlin, 10099 Berlin, Germany

###### Abstract

Extreme multistability (EM) is characterized by the emergence of infinitely
many coexisting attractors or continuous families of stable states in
dynamical systems. EM implies complex and hardly predictable asymptotic
dynamical behavior. We analyse a model for pendulum clocks coupled by springs
and suspended on an oscillating base, and show how EM can be induced in this
system by a specifically designed coupling. First, we uncover that symmetric
coupling can increase the dynamical complexity. In particular, the coexistence
of multiple isolated attractors and continuous families of stable periodic
states is generated in a symmetric cross-coupling scheme of four pendulums.
These coexisting infinitely many states are characterized by different levels
of phase synchronization between the pendulums, including anti-phase and in-
phase states. Some of the states are characterized by splitting of the
pendulums into groups with silent sub-threshold and oscillating behavior,
respectively. The analysis of the basins of attraction further reveals the
complex dependence of EM on initial conditions.

††preprint: AIP/123-QED

> The coexistence of several asymptotic stable states for a dynamical system
> with fixed parameter values is called multistability. This phenomenon has
> been identified in diverse fields of science both experimentally and
> theoretically. Which asymptotic state the system will converge to is
> determined solely by its initial state. When the number of stable states is
> infinite, extreme multistability (EM) becomes a dominant feature.
> Understanding EM and its control is an important issue, because systems with
> EM offer even greater flexibility than those with finite multistability when
> switching from one stable state to another. We give an example of EM in a
> coupled pendulums model that takes into account an escapement mechanism as
> well as local and global couplings. We have paid a particular attention to
> the coupling structure that leads to the emergence of EM.

##  I Introduction

Complex networks have largely enriched our understanding of a variety of
complex dynamical systems in many fields, such as biology, ecology,
climatology, sociology, and others Newman (2018); Yanchuk  _et al._ (2021). By
modeling real-world systems as networks in which collections of dynamic nodes
are connected by static or adaptive edges, one can study collective behaviors
both analytically and numerically Blekhman (1988); Arenas  _et al._ (2008);
Newman (2003); Gross and Blasius (2008); Yanchuk  _et al._ (2022).

Synchronization is a ubiquitous dynamical phenomenon that has been observed in
many natural and engineering systems Arenas  _et al._ (2008); Pikovsky  _et
al._ (2001); Osipov  _et al._ (2007); Boda  _et al._ (2013); Chhabria  _et
al._ (2018); Ikeguchi and Shimada (2019). Different types of synchronous
patterns have been identified involving complete synchronization Pecora and
Carroll (1998) (oscillators’ states become asymptotically the same with time),
cluster synchronization (a network splits into groups of synchronous elements)
Dahms  _et al._ (2012); Lücken and Yanchuk (2012), special types of spatial
coexistence of coherent and incoherent states Abrams and Strogatz (2004);
Omelchenko  _et al._ (2011); Kasatkin  _et al._ (2017), and many others.
Various patterns have been found in experimental contexts, such as
optoelectronic networks Soriano  _et al._ (2013), chemical networks Tinsley
_et al._ (2012), neural networks Hammond  _et al._ (2007), ecological Blasius
_et al._ (1999), and climate systems Tsonis and Roebber (2004).

Apart from synchronization, multistability – the coexistence of several
asymptotic stable states (attractors) for a given set of parameters – is
another intriguing phenomenon which has been studied for decades in modern
nonlinear science Feudel (2008); Pisarchik and Feudel (2014); Dudkowski  _et
al._ (2022). The final state of a system with multistability depends crucially
on initial conditions. Multistability has also been observed in many areas of
science, such as nonlinear optics Arecchi (1991), neuroscience Foss  _et al._
(1996), climate dynamics Power and Kleeman (1993), laser physics Masoller
(2002), electronic oscillators Borresen and Lynch (2002), and in different
classes of systems, such as weakly dissipative systems Feudel  _et al._
(1996), systems with time delays Yanchuk and Wolfrum (2010); Balanov  _et al._
(2005), and coupled systems Feudel  _et al._ (1998).

Understanding the emergence of co-existing attractors is an important issue,
and controlling multistability is an even more difficult task. When the number
of co-existing attractors increases infinitely, EM emerges. In coupled
systems, the presence of EM has been found to be closely related to partial
synchrony Hens  _et al._ (2015). By designing a specific coupling scheme to
achieve partial synchrony, one can obtain infinitely many coexisting stable
states Hens  _et al._ (2015); Sun  _et al._ (1999); Ngonghala  _et al._
(2011); Hens  _et al._ (2012); Pal  _et al._ (2014). Apart from the
conservative cases, a common reason for the occurrence of EM in networks is
time-reversibility, a special type of spatio-temporal symmetry Politi  _et
al._ (1986); Pikovsky and Rosenau (2006); Lamb and Roberts (1998); Ashwin  _et
al._ (2016); Burylko  _et al._ (2018).

Despite the extensive literature on multistable dynamical systems, the
emergence of multistability or EM in networked dynamical systems remains a
challenging problem due to a large number of possible routes to EM, some of
which have yet to be discovered. Analytical and numerical challenges arise
from the diversity of coupling topologies and the complexity of individual
models.

In this work, we address the multistability problem in a mathematical model of
coupled clocks suspended on a rotating disc and additionally coupled with
springs. The interaction of the clocks with the disc provides the global
coupling among all clocks and therefore influences their behavior, similar to
the interaction of the pedestrians with the bridge in the famous effect of
crowd synchrony on the Millennium Bridge Strogatz  _et al._ (2005). Such a
global scheme has also proved useful in uncovering complex transient states
Dudkowski  _et al._ (2020a). The oscillating clocks are also locally coupled
via springs. In Ref. Dudkowski  _et al._ , 2020b, a similar system of three
coupled clocks was studied.

The following main results are obtained in this work:

  * •

We generalize the system of three coupled clocks Dudkowski  _et al._ (2020b)
into a network-coupled scenario allowing arbitrary coupling configurations.

  * •

We investigate how different coupling topologies affect the multistability in
systems of three and four coupled clocks. We observe that more symmetric
coupling topologies can lead to more complex dynamics with higher
multistability.

A particularly reach appears to be the “cross-coupling” structure with
“diagonal” spring couplings in the system of four coupled clocks. In such a
case, we observe EM that combines continuous family of stable attractors with
different phase relations between the clocks. We provide an analytical and
numerical description of this new phenomenon.

  * •

Furthermore, we discuss how the discontinuity of the escapement mechanism
affects the multistability in the system. The clocks within certain coupled
groups (clusters) remain either silent or oscillating and in-phase
synchronized due to the switching of the escapement mechanism. This leads to
three qualitatively different discontinuity-induced types of attractors.

Table 1: Parameters for the $N$-pendulums of system (1). Parameter | Value |  Definition
---|---|---
$B_{0}$ | 5.115 [kgm2] |  Support’s moment of inertia
$k_{\varphi}$ | 17.75 [N/m] |  Stiffness coefficient of springs
$k_{\theta}$ | 34 [Nm] |  Stiffness of the spring connecting the base and the unmoving support
$c_{\theta}$ | $\Delta=ln(2)$ |  Damping of the supporting base
$c_{\varphi}$ | 0.01 [Nms] |  Damping of the damper connecting the base and the unmoving support
$m$ | 1 [kg] |  Mass of each pendulum clock
$l$ | 0.24849 [m] |  Lengths of the pendulums
$g$ | 9.81 [m/s2] |  Gravity acceleration
$d$ | 1 [m] |  Distance between $O$ and $S_{i}$ $(i=1,2,...,N)$
$\alpha_{i}=\sphericalangle({Ox},~{}\overline{OS_{i}})$ | $(\frac{360}{N})^{\circ}=\frac{2\pi}{N}$ |  $Ox$ denotes the positive $x$ half axis
$M$ | 0.075 [Nm] |  Fixed external momentum
$\varepsilon_{0}$ | $5^{\circ}=5*\pi/180$ |  Escapement threshold

##  II Model and measures

###  II.1 General model

We first present a mathematical model of the $N$ coupled pendulum clocks
suspended on a rotating disc, see Fig. 1(A). The rotating disc provides a
global coupling, while the springs allow for arbitrary local coupling
structure. Our model is a generalization of the system of three pendulums from
Ref. Dudkowski  _et al._ , 2020b.

Figure 1:  (A) The scheme of $N$ coupled identical pendulum clocks (shown as
circles with different colors) suspended at evenly distributed black points
$S_{i}$, $i=1,2,...,N$, on an oscillating supporting base (the $xy$ plane).
The local coupling is realized using springs between the clocks. (B)–(D): For
$N=3$, three types of coupling structures of springs include all-to-all,
asymmetric, and symmetric topologies. (E)–(G): For $N=4$, three types of
coupling structures of springs include all-to-all, asymmetric, and symmetric
topologies.

The supporting base is placed at the origin $O$ of the $xy$ plane and it can
oscillate freely around the axis perpendicular to the plane of Fig. 1; the
angular deviation of the base is $\theta$. The properties of the base are
described by the moment of inertia $B_{0}$ [kgm2], the stiffness of the spring
connecting the base and the static support $k_{\theta}$ [Nm], and the damping
$c_{\theta}$ [Nms]. Identical pendulums (marked colored filled circles) are
suspended at evenly distributed black points $S_{i}$, $i=1,2,3,...,N$, i.e.,
the angles between ${OS}_{i}$ and ${OS}_{i+1}$ (index $i$ is considered mod
$N$) are ${{2\pi}/{N}}$. The angles
$\alpha_{i}=\sphericalangle({Ox},~{}\overline{OS}_{i})$ characterize the
angular position of the suspension points $S_{i}$, where $Ox$ is the positive
$x$ half-axis. The parameter $d=|\overline{OS}_{i}|$ $(i=1,2,...,n)$ is the
distance between the origin $O$ and each suspension point $S_{i}$. Each
pendulum is described by the angle displacement $\varphi_{i}$, the mass $m$
[kg], the length $l$ [m], and the damping coefficient $c_{\varphi}$ [Nms]. The
stiffness coefficients of the springs are $k_{\varphi}$ [N/m]. The description
and the values for all parameters are summarized in Table 1.

The equations of motion of the $N$ coupled pendulums is given by the following
system:

$\begin{array}[]{ll}\displaystyle(B_{0}+nmr^{2})\ddot{\theta}+k_{\theta}\theta+c_{\theta}\dot{\theta}+\sum\limits_{i=1}^{n}mr\\{l[\ddot{\varphi}_{i}\sin(\varphi_{i}-\theta-\alpha_{i})+\dot{\varphi}^{2}_{i}\cos(\varphi_{i}-\theta-\alpha_{i})]+g\cos(\alpha_{i}+\theta)\\}+\Delta
V_{\theta}=0,\\\
ml^{2}\ddot{\varphi}_{i}+mgl\sin{\varphi_{i}}+c_{\varphi}{\dot{\varphi}_{i}}+mrl[\ddot{\theta}\sin(\varphi_{i}-\theta-\alpha_{i})-\dot{\theta}^{2}\cos(\varphi_{i}-\theta-\alpha_{i})]+\Delta
V_{\varphi_{i}}=M_{E_{i}},\\\ \end{array}$ (1)

where $i=1,2,...,N$.

The build-in escapement mechanism produces the moment of force, which is
modeled by the discontinuous functions $M_{E_{i}}$, $i=1,2,3,...,N$
Czołczyński  _et al._ (2011); Kapitaniak  _et al._ (2012). These functions
depend not only on the displacement $\varphi_{i}(t)$, but also on the position
of the $i$-th mechanism’s cogwheel versus the mechanism’s pallet
$\sigma_{i}(t)$:

$M_{E_{i}}=\begin{cases}M&:\sigma_{i}=1\wedge
0<\varphi_{i}<\varepsilon_{0},\\\
-M&:\sigma_{i}=2\wedge-\varepsilon_{0}<\varphi_{i}<0,\\\
0&:\mathrm{otherwise.}\end{cases}$ (2)

Here $M=0.075$ [Nm] represents the value of the external momentum, while
$\varepsilon_{0}=5.0^{\circ}$ denotes the escapement’s threshold (the
mechanism turns off as the pendulum exceeds this threshold). In fact,
$\sigma_{i}(t)$ become additional discrete-valued variables in the system that
are influencing the system’s dynamics via the terms $M_{E_{i}}$ and which are
changing discontinuously according to the following rules:

  * (I)

When a pendulum $\varphi_{i}$ crosses the escapement threshold at some time
moment $t^{*}$: $\varphi_{i}(t_{*})=\varepsilon_{0}$ with increasing
$\varphi_{i}$, i.e., $\dot{\varphi}_{i}(t_{*})>0$, the variable
$\sigma_{i}(t)$ is set to 2 for all $t\in[t^{*},t_{\mathrm{e}})$, where
$t_{\mathrm{e}}$ is the time of a next event.

  * (II)

When the pendulum $\varphi_{i}$ crosses the escapement threshold
$\varphi_{i}(t_{*})=-\varepsilon_{0}$ with decreasing $\varphi_{i}$, i.e.,
$\dot{\varphi}_{i}(t_{*})<0$, the variable $\sigma_{i}(t)$ is set to 1 for all
$t\in[t^{*},~{}t_{\mathrm{e}})$, where $t_{\mathrm{e}}$ is the time moment of
a next crossing event.

In this way, the variables $\sigma_{i}(t)$ are piece-wise constant with the
possible discrete values 1 or 2, which change discontinuously when either
event (I) or (II) occurs.

The terms $\Delta V_{\varphi_{i}}$ and $\Delta V_{\theta}$ in model (1)
describe the moments of forces from the coupling springs. These two terms can
be written explicitly using the following terms: $s_{ij}$, the constant
distance between the $i$-th and $j$-th clocks when the system stays still, and
$\hat{s}_{ij}(t)$, the time-dependent distance between the $i$-th and $j$-th
clocks for the moving system:

$\begin{array}[]{ll}\displaystyle
s_{ij}=r\sqrt{2(1-\cos(\alpha_{i}-\alpha_{j}))},\\\
\hat{s}_{ij}=\sqrt{s_{ij}^{2}+2l^{2}(1-\cos(\varphi_{i}-\varphi_{j}))+8lr\sin\left(\frac{\varphi_{i}-\varphi_{j}}{2}\right)\sin\left(\frac{\alpha_{i}-\alpha_{j}}{2}\right)\sin\left(\frac{\varphi_{i}+\varphi_{j}-\alpha_{i}-\alpha_{j}}{2}-\theta\right)},\\\
\Delta
V_{\theta}=2lrk_{\varphi}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}a_{ij}\left(1-\frac{s_{ij}}{\hat{s}_{ij}}\right)\sin\left(\frac{\varphi_{i}-\varphi_{j}}{2}\right)\sin\left(\frac{\alpha_{j}-\alpha_{i}}{2}\right)\cos\left(\frac{\varphi_{i}+\varphi_{j}-\alpha_{i}-\alpha_{j}}{2}-\theta\right),\\\
\Delta
V_{\varphi_{i}}=\sum\limits_{i=1}^{n}a_{ij}k_{\varphi}l\left(1-\frac{s_{ij}}{\hat{s}_{ij}}\right)\left[l\sin(\varphi_{i}-\varphi_{j})+2r\sin\left(\frac{\alpha_{i}-\alpha_{j}}{2}\right)\sin\left(\varphi_{i}-\frac{\alpha_{i}+\alpha_{j}}{2}-\theta\right)\right],\end{array}$
(3)

where $(a_{ij})$ is the coupling matrix via the springs, i.e., $a_{ij}=1$ if
the pendulum $i$ is connected with the pendulum $j$ via a spring and
$a_{ij}=0$ otherwise. $a_{ii}=0$ since there are no self-loops.

Figure 2:  Distributions of the order parameter $r$ for different coupling
topologies. (A)-(C) show for 3-coupled pendulums ($N=3$), the identified
multistabilities for all-to-all, asymmetric, and symmetric coupling structures
(see right-upper corners), respectively. For each coupling topology, 5000
order parameters are obtained from 5000 simulations using random initial
conditions. Parameters for $N=3$ are fixed as in Table 1, in particular,
$\alpha_{1}=\frac{\pi}{2},$ $\alpha_{2}=\frac{7\pi}{12}$, and
$\alpha_{3}=\frac{11\pi}{12}$. Each vertical bar corresponds to a potential
attractor. Similarly, (D)-(F) are multistabilities for 4-coupled pendulums
($N=4$), based on 1000 order parameters obtained from 1000 simulations.
Parameters for $N=4$ are also fixed as in Table 1, in particular,
$\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and
$\alpha_{4}=2\pi$. For $N=3$, compared with (A) and (C), the asymmetric
coupling structure in (B) decreases the dynamical complexity; while for $N=4$,
compared with (E), (F) shows that the symmetric coupling can increase the
dynamical complexity due to the emergence of EM.

The influence of different coupling structures on collective dynamics has not
been systematically reported for this model. In the remaining part of this
paper, we consider the cases $N=3$ (Figs. 1(B)-(D)) and $N=4$ (Figs.
1(E)-(G)). In particular, we focus on the following questions:

  * (1)

How do different coupling topologies alter synchronization states and their
basins of attraction?

  * (2)

Is there extreme multistability in the coupled pendulum model? If so, what is
its origin?

  * (3)

What are the effects of the discontinuity of the escapement mechanism on the
dynamics and multistability?

We fix the parameters as in Table 1. For $N=3$, the three considered types of
coupling structures of springs include all-to-all (Fig. 1(B)), asymmetric
(Fig. 1(C)), and symmetric (Fig. 1(D); mirror symmetry with respect to the
vertical axis) topologies. Identical pendulums are suspended at evenly
distributed points $S_{i}$, $i=1,2,3$ with
$\sphericalangle(\overline{OS}_{i},~{}\overline{OS}_{i+1})=120^{\circ}$. For
$N=4$, we also consider three types of coupling structures of springs: all-to-
all (Fig. 1(D)), asymmetric (Fig. 1(E)), and symmetric (Fig. 1(F)) topologies.
Here also the identical pendulums are suspended at evenly distributed points
$S_{i}$, $i=1,2,3,4$ with
$\sphericalangle(\overline{OS}_{i},~{}\overline{OS}_{i+1})=90^{\circ}$.

We use Monte Carlo sampling and two classical measures (order parameter and
mean frequencies) for the numerical quantification of synchronization states
and the analysis of their basins of attraction.

###  II.2 Measures

Order parameter. We visualize the dynamics of the synchronization transitions
with the Kuramoto order parameter:

$R(t)=\frac{1}{N}\sum\limits_{i=1}^{N}e^{i\varphi_{i}(t)},$ (4)

where $N$ is the number of oscillators. When $|R(t)|=1$ ($|R(t)|\approx 0$),
oscillators are in the complete synchronization (disordered) state. The degree
of synchronization in numerical simulations is quantified using the averaged
value of the order parameter:

$r=\frac{1}{T_{\mathrm{av}}}\int_{T_{\mathrm{tr}}}^{T_{\mathrm{tr}}+T_{\mathrm{av}}}|R(t)|dt.$
(5)

over the time interval $T_{\mathrm{av}}=50$ after a sufficiently long
transient time $T_{\mathrm{tr}}$.

Mean frequency. The mean oscillation frequency of a pendulum is given as:

$\langle\omega_{i}\rangle=\frac{2\pi n_{i}}{T_{\mathrm{av}}},$ (6)

where the same time interval of $T_{\mathrm{av}}=50$ is applied and $n_{i}$
represents the number of complete oscillations of $i$th clock within this
interval. The number of complete oscillations can be computed using the number
of crossings of the Poincare map $\varphi_{i}=0$ or
$\varphi_{i}=\varepsilon_{0}$. The mean frequency is calculated using the last
$T_{\mathrm{av}}=50$ time units after a sufficiently long transient time
$T_{\mathrm{tr}}$.

##  III Collective dynamics for different coupling topologies

We first conduct various simulations of the system for random initial
conditions. More specifically, we choose the following initial conditions
$[\theta^{0}=0.01,\phi_{1}^{0},\dots,\phi_{n}^{0},\dot{\theta^{0}}=0,\dot{\phi_{1}^{0}}=0,\dots,\dot{\phi_{n}^{0}}=0]$,
where $\phi_{1}^{0},\dots,\phi_{n}^{0}$, are chosen randomly from the interval
$[-\pi,\pi)$. For the case $N=3$, simulations with 5,000 different initial
conditions are performed with the integration time 15,000, and the last 50
time units are used for the calculation of the order parameter and the mean
frequency. For $N=4$, we perform 1,000 runs with the integration interval
10,000, and the transient 9,950. The $r$ and $\langle\omega_{i}\rangle$ from
Eqs. (5) and (6), respectively, are used to estimate the synchronization state
(attractor) in each simulation. We found that further increase of the number
of runs and the integration interval does not affect the obtained results
qualitatively.

Figure 2 shows the distribution of the order parameter $r$ for different
initial conditions. This distribution reveals the possible number of different
attractors. Figures 2(A)-(C) correspond to the coupling structures of three
clocks in the Figs. 1(B)-(D), respectively. Here we see finitely many isolated
lines indicating a relatively small number of possible synchronization states.
Interestingly, the case of asymmetric coupling in Fig. 1(C) exhibits lower
dynamical complexity as that shown in Fig. 2(B), since only one line of the
distribution of $r$ is achieved for all initial conditions. Figures 2(A) and
2(C) imply finite multistability with different possible asymptotic values of
$r$.

For four coupled clocks (Figs. 2(D)-(F)), we also uncover that different
coupling topologies lead to diverse dynamical complexities. Specifically,
three different lines of $r$ are observed in Figs. 2(D) and 2(E). More
importantly, if the structure of the coupled clocks maintains the symmetry as
Fig. 1(G), the distribution of the asymptotic order parameters in Fig. 2(F) is
no longer discrete, but contains continuous parts. Such a distribution
indicates higher complexity and even EM. In order to characterize deeper the
emergence of EM, we focus on analytical and numerical explanations of this
phenomenon in the following sections.

##  IV Extreme multistability

We recall that EM is potentially observed for the cross-coupling structure of
four coupled clocks (Fig. 1(G)), where the distribution of asymptotic order
parameters seems to be continuous (Fig. 2(F)). In this scheme, the opposite
pendulums are connected by springs. The corresponding coupling matrix has four
nonzero entries $a_{13}=a_{31}=a_{24}=a_{42}=1$, and the angle position
parameters are $\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$,
$\alpha_{3}=\frac{3\pi}{2}$ and $\alpha_{4}=2\pi$. The system (1) becomes:

$\displaystyle\begin{array}[]{ll}\displaystyle(B_{0}+4mr^{2})\ddot{\theta}+k_{\theta}\theta+c_{\theta}\dot{\theta}+mrl[-\ddot{\varphi}_{1}\cos(\theta-\varphi_{1})-\dot{\varphi}^{2}_{1}\sin(\theta-\varphi_{1})]+mrl[\ddot{\varphi}_{2}\sin(\theta-\varphi_{2})-\dot{\varphi}^{2}_{2}\cos(\theta-\varphi_{2})]\\\
+mrl[\ddot{\varphi}_{3}\cos(\theta-\varphi_{3})+\dot{\varphi}^{2}_{3}\sin(\theta-\varphi_{3})]+mrl[-\ddot{\varphi}_{4}\sin(\theta-\varphi_{4})+\dot{\varphi}^{2}_{4}\cos(\theta-\varphi_{4})]+\Delta
V_{\theta}=0,\\\
ml^{2}\ddot{\varphi}_{1}+mgl\sin{\varphi_{1}}+c_{\varphi}{\dot{\varphi}_{1}}+mrl[-\ddot{\theta}\cos(\theta-\varphi_{1})+\dot{\theta}^{2}\sin(\theta-\varphi_{1})]+\Delta
V_{\varphi_{1}}=M_{E_{1}},\\\
ml^{2}\ddot{\varphi}_{2}+mgl\sin{\varphi_{2}}+c_{\varphi}{\dot{\varphi}_{2}}+mrl[\ddot{\theta}\sin(\theta-\varphi_{2})+\dot{\theta}^{2}\cos(\theta-\varphi_{2})]+\Delta
V_{\varphi_{2}}=M_{E_{2}},\\\
ml^{2}\ddot{\varphi}_{3}+mgl\sin{\varphi_{3}}+c_{\varphi}{\dot{\varphi}_{3}}+mrl[\ddot{\theta}\cos(\theta-\varphi_{3})-\dot{\theta}^{2}\sin(\theta-\varphi_{3})]+\Delta
V_{\varphi_{3}}=M_{E_{3}},\\\
ml^{2}\ddot{\varphi}_{4}+mgl\sin{\varphi_{4}}+c_{\varphi}{\dot{\varphi}_{4}}+mrl[-\ddot{\theta}\sin(\theta-\varphi_{4})-\dot{\theta}^{2}\cos(\theta-\varphi_{4})]+\Delta
V_{\varphi_{4}}=M_{E_{4}},\\\ \end{array}$ (13)

where

$\begin{array}[]{ll}\displaystyle\Delta
V_{\theta}=4lrk_{\varphi}\left[\left(1-\frac{s_{24}}{\hat{s}_{24}}\right)\sin\left(\frac{\varphi_{2}-\varphi_{4}}{2}\right)\sin\left(\theta-\frac{\varphi_{2}+\varphi_{4}}{2}\right)-\left(1-\frac{s_{13}}{\hat{s}_{13}}\right)\sin\left(\frac{\varphi_{1}-\varphi_{3}}{2}\right)\cos\left(\theta-\frac{\varphi_{1}+\varphi_{3}}{2}\right)\right],\\\
\Delta
V_{\varphi_{1}}=k_{\varphi}l\left(1-\frac{2r}{\hat{s}_{13}}\right)[l\sin(\varphi_{1}-\varphi_{3})-2r\sin(\theta-\varphi_{1})],\\\
\Delta
V_{\varphi_{2}}=k_{\varphi}l\left(1-\frac{2r}{\hat{s}_{24}}\right)[l\sin(\varphi_{2}-\varphi_{4})-2r\cos(\theta-\varphi_{2})],\\\
\Delta
V_{\varphi_{3}}=k_{\varphi}l\left(1-\frac{2r}{\hat{s}_{31}}\right)[l\sin(\varphi_{3}-\varphi_{1})+2r\sin(\theta-\varphi_{3})],\\\
\Delta
V_{\varphi_{4}}=k_{\varphi}l\left(1-\frac{2r}{\hat{s}_{24}}\right)[l\sin(\varphi_{4}-\varphi_{2})+2r\cos(\theta-\varphi_{4})].\end{array}$
(14)

The distances $s_{ij}$ and $\hat{s}_{ij}$ are:

$\begin{array}[]{ll}\displaystyle s_{13}=s_{31}=s_{24}=s_{42}=2r,\\\
\hat{s}_{13}=\hat{s}_{31}=\sqrt{4r^{2}+2r^{2}(1-\cos(\varphi_{1}-\varphi_{3}))-8lr\sin\left(\frac{\varphi_{1}-\varphi_{3}}{2}\right)\sin\left(\theta-\frac{\varphi_{1}+\varphi_{3}}{2}\right)},\\\
\hat{s}_{24}=\hat{s}_{42}=\sqrt{4r^{2}+2r^{2}(1-\cos(\varphi_{2}-\varphi_{4}))-8lr\sin\left(\frac{\varphi_{2}-\varphi_{4}}{2}\right)\cos\left(\theta-\frac{\varphi_{2}+\varphi_{4}}{2}\right)}.\\\
\end{array}$ (15)

As we will see later, the regime of EM is characterized by the emergence of
two frequency synchronized clusters each containing two clocks. The following
phase relations are observed for the synchronized clusters: “in-phase-in-
phase" (II), “in-phase-anti-phase" (IA), “anti-phase-in-phase" (AI) and “anti-
phase-anti-phase" (AA). The exact meaning of these relations are given in
Table 2. For example, IA means that the clocks in the first cluster are in-
phase and anti-phase in the second cluster. Additionally, due to the
discontinuity induced by the escapement mechanism, the mixed states are
observed, when one or both of the clusters are not oscillating. This is
possible due to the fact that the clocks do not cross periodically the
escapement threshold and, hence, do not gain energy. The following clusters
are observed: "silent-in-phase" (SI), "in-phase-silent" (IS), and "silent-
silent" (SS).

Figure 3:  Distribution of the order parameter $r$ for 4 coupled clocks (13) with the coupling topology as in Fig. 1(G). (A) is same as Fig. 2(F), obtained by Monte Carlo sampling with 1,000 random trials. (B)-(E) represent the parts of the distribution (extracted from (A)) that correspond to specific cluster states: (B) counts only the order parameters for the trials ending in II (in-phase-in-phase) configuration, (C) stands for IA, (D) for AI, and (E) for AA, see Table 2 explaining the cluster states. The main observation is that only case (B) is related to the emergence of EM. Parameters are fixed as in Table 1 with $\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and $\alpha_{4}=2\pi$ for the 4-coupled clocks ($N=4$).  Table 2: Different phase-clusters possibilities for 4-coupled clocks (13) with cross-coupling structure (see Fig. 1(G)). Symbol | Description
---|---
II | In-phase-in-phase ($\varphi_{1}(t)=\varphi_{3}(t),\ \varphi_{2}(t)=\varphi_{4}(t)$)
IA | In-phase-anti-phase ($\varphi_{1}(t)=\varphi_{3}(t),\ \varphi_{2}(t)=-\varphi_{4}(t)$)
AI | Anti-phase-in-phase ($\varphi_{1}(t)=-\varphi_{3}(t),\ \varphi_{2}(t)=\varphi_{4}(t)$)
AA | Anti-phase-anti-phase ($\varphi_{1}(t)=-\varphi_{3}(t),\ \varphi_{2}(t)=-\varphi_{4}(t)$)
SI | Silent-in-phase ($\varphi_{1}(t)=\varphi_{3}(t)=0,\ \varphi_{2}(t)=\varphi_{4}(t)$)
IS | In-phase-silent ($\varphi_{1}(t)=\varphi_{3}(t),\ \varphi_{2}(t)=\varphi_{4}(t)=0$)
SS | Silent-silent ($\varphi_{1}(t)=\varphi_{3}(t)=0,\ \varphi_{2}(t)=\varphi_{4}(t)=0$)

In Figure 3, we split the probability distribution of the order parameter $r$
accordingly to the cluster states observed. Specifically, Fig. 3(A) gives the
whole distributions, same as in Fig. 2(F)). Figure 3(B) exacts from Fig. 3(A)
only the order parameters that correspond to II phase clusters, Fig. 3(C) to
IA, Fig. 3(D) to AI, and Fig. 3(E) to AA. Only II clusters exhibit a
continuous distribution of $r$, thus suggesting that EM appears due to such
type of clusters. In contrast, Figs. 3(C)-(E) shows only a finite number of
lines of $r$. In the following sections, we provide additional analytical and
numerical evidences that confirm our observation and explain the phenomenon of
EM.

###  IV.1 Family of stable cluster states

####  IV.1.1 Theoretical analysis of EM

As Fig. 3(B) indicates, coexistence of infinitely many stable states can be
related to the emergence of II, the in-phase-in-phase clusters. To study the
existence of such clusters, we show that the following subspace of II
solutions:

$\varphi_{1}(t)=\varphi_{3}(t)=\psi_{1}(t),\quad\varphi_{2}(t)=\varphi_{4}(t)=\psi_{2}(t),\quad\theta(t)=0.$
(16)

is invariant with respect to the solutions of system (13). Indeed,
substituting $\varphi_{1}(t)=\varphi_{3}(t)=\psi_{1}(t)$ and
$\varphi_{2}(t)=\varphi_{4}(t)=\psi_{2}(t)$ into (13), we obtain $\Delta
V_{\theta}=0$ and $\Delta V_{\varphi_{i}}=0$, and the equations for the new
variables $\theta$, $\psi_{1}$, and $\psi_{2}$ read:

$\displaystyle(B_{0}+4mr^{2})\ddot{\theta}+k_{\theta}\theta+c_{\theta}\dot{\theta}=0,$
(17a) $\displaystyle
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}+mrl[-\ddot{\theta}\cos(\theta-\psi_{1})+\dot{\theta}^{2}\sin(\theta-\psi_{1})]=M_{E_{1}},$
(17b) $\displaystyle
ml^{2}\ddot{\psi}_{2}+mgl\sin{\psi_{2}}+c_{\varphi}{\dot{\psi}_{2}}+mrl[\ddot{\theta}\sin(\theta-\psi_{2})+\dot{\theta}^{2}\cos(\theta-\psi_{2})]=M_{E_{2}},$
(17c) $\displaystyle
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}+mrl[\ddot{\theta}\cos(\theta-\psi_{1})-\dot{\theta}^{2}\sin(\theta-\psi_{1})]=M_{E_{1}},$
(17d) $\displaystyle
ml^{2}\ddot{\psi}_{2}+mgl\sin{\psi_{2}}+c_{\varphi}{\dot{\bar{\varphi}}_{2}}+mrl[-\ddot{\theta}\sin(\theta-\psi_{2})-\dot{\theta}^{2}\cos(\theta-\psi_{2})]=M_{E_{2}},$
(17e)

with

$M_{E_{i}}=\begin{cases}M&:\sigma_{i}=1\wedge 0<\psi_{i}<\varepsilon_{0}\\\
-M&:\sigma_{i}=2\wedge-\varepsilon_{0}<\psi_{i}<0\\\
0&:\mathrm{otherwise}\end{cases}$ (18)

where $i=1,2$, $M=0.075$ [Nm], and $\varepsilon=5.0^{\circ}$.

Now, by setting $\theta=0$, we observe that Eq. (17a) is satisfied, and the
dynamical equations for $\psi_{1}$ (Eqs. (17b) and (17d)) and for $\psi_{2}$
(Eqs. (17c) and (17e)) become the same:

$\displaystyle
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}=M_{E_{1}},$
(19a) $\displaystyle
ml^{2}\ddot{\psi}_{2}+mgl\sin{\psi_{2}}+c_{\varphi}{\dot{\psi}_{2}}=M_{E_{2}}.$
(19b)

Thus, we have proven the following:

######  Proposition 1 (II-cluster subspace)

The subspace (16) of the cluster II solutions is invariant with respect to the
solutions (flow) of system (13). This subspace is 4-dimensional
($S^{2}\times\mathrm{R}^{2}$), and the flow on this subspace is given by Eqs.
(19a) and (19b) which describe the relative motion of the clusters.

Another important observation is that the dynamical equations (19a) for
$\psi_{1}$ and (19b) $\psi_{2}$ in the II-cluster subspace are (i) the same,
and (ii) uncoupled from each other, i.e., the equation for $\psi_{1}$ is
independent on $\psi_{2}$ and vise-versa. The latter property leads to the
coexistence of infinitely many asymptotic states (and to EM eventually). We
formulate the corresponding result as a proposition.

######  Proposition 2 (EM of II-cluster states)

Assume that $c_{\varphi},m,l,g,k_{\theta},c_{\theta}$ and $B_{0}$ are positive
parameters. Assume also that system (19a) (or, equivalently, (19b)) possesses
a stable nontrivial asymptotic state (attractor), and $\psi^{*}(t)$ is a
solution on this attractor (i.e., a single clock has a stable oscillatory
state). Then the system on the II-cluster subspace (19a)–(19b) possesses the
following asymptotic states:
– II clusters:

$\varphi_{1}=\varphi_{3}=\psi_{1}=\psi^{*}(t),\quad\varphi_{2}=\varphi_{4}=\psi_{2}=\psi^{*}(t+\gamma).$
(20)

– SI clusters:

$\varphi_{1}=\varphi_{3}=\psi_{1}=0,\quad\varphi_{2}=\varphi_{4}=\psi_{2}=\psi^{*}(t).$
(21)

– SS clusters:

$\varphi_{1}=\varphi_{3}=\psi_{1}=0,\quad\varphi_{2}=\varphi_{4}=\psi_{2}=0.$
(22)

– IS clusters:

$\varphi_{1}=\varphi_{3}=\psi_{1}=\psi^{*}(t),\quad\varphi_{2}=\varphi_{4}=\psi_{2}=0.$
(23)

where $\gamma$ is an arbitrary real constant describing a phase shift between
the clusters. Moreover, if $\gamma^{*}(t)$ is an orbitally asymptotically
stable limit cycle (stable periodic oscillations of the clock), then the
states (20) build a stable invariant torus foliated by (infinitely many)
periodic solutions of the form (20).

The main message of Proposition 2 is that under "normal conditions" when the
single clock oscillates periodically, the coupled system can have stable
cluster II oscillations with an arbitrary phase-shift between the clusters. If
the phase-shift is zero, the order parameter $r$ is highest and equal 1, while
it can achieve a continuous range of smaller values depending on the phase
shift. The coexistence of such states leads to EM. Note that the states SI,
SS, and IS do not lead to EM, but correspond to isolated attractors in the
coupled system.

Proof of Proposition 2. Let us first mention that the equilibrium
$\psi_{1}=\psi_{2}=0$ is asymptotically stable, and it corresponds to the
stability of a silent state of a pendulum with damping and without external
energy inflow. Therefore, the two independent systems (19a) and (19b) can
reach both attractors: zero equilibrium and non-trivial attractor
corresponding to $\psi^{*}(t)$, depending on initial conditions. Moreover, an
arbitrary phase shift $\psi^{*}(t+\gamma)$ is clearly also possible and belong
to the same nontrivial attractor. This provides the existence of the states
(20)–(23).

The invariant torus from the proposition corresponds to the direct product of
the limit cycles $\mathcal{C}\times\mathcal{C}$, where
$\mathcal{C}=\\{(\psi,~{}\dot{\psi})\in(S^{1}\times\mathbb{R}):\psi=\gamma^{*}(t),~{}t\in\mathbb{R}\\}$.
The stability of this torus follows from the orbital stability of the limit
cycle in each of the subsystem and the properties of the cross-product of the
uncoupled system (19a)–(19b). End of proof.

Figure 4:  Three exemplary dynamical patterns for the regime of EM (see Fig.
3(B)) for 4-coupled pendulums with cross-coupling structure (13). (A) Almost
complete synchronization with a small phase shift between the clusters of
clocks. (B) Phase synchronization with an intermediate phase shift between the
clusters. (C) Anti-phase synchronization between the clusters. Details of
column figures are as follows: a: the time series for variables $\theta$,
$\varphi_{1}$, $\varphi_{2}$, $\varphi_{3}$, $\varphi_{4}$, respectively, b:
the phase-time plots of the pendula, c: the mean frequencies of the clocks, d:
Projections on two phase variables form different clusters (orange lines) and
Poincaré maps (black dots). Parameters are fixed as in Table 1 with
$\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and
$\alpha_{4}=2\pi$ for the 4-coupled clocks ($N=4$).  Figure 5:  Dynamical
patterns of three types of isolated attractors with large-amplitude
oscillations for 4 coupled pendulums with cross-coupling structure (13).
(A)-(C): IA, AI, and AA synchronization patterns (see Table 2 for the
explanation of the abbreviations). IA, AI, and AA correspond to the order
parameters from Figs. 3(C), (D), and (E), respectively. Information on the
columns: a: the time series for variables $\theta$, $\varphi_{1}$,
$\varphi_{2}$, $\varphi_{3}$, $\varphi_{4}$, respectively; b: phase-time plots
of the pendulum angles; c mean frequencies; d: projection of the solution on
the phase variables from different clusters (orange lines) and Poincaré maps
(black dots). Parameters are fixed as in Table 1 with
$\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and
$\alpha_{4}=2\pi$ for the 4-coupled clocks ($N=4$).

####  IV.1.2 Numerical study of EM

Figure 4 shows three examples of different stable II synchronization patterns
from the continuous family of solutions by Eq. (20). Figure 4 provides (A)
almost in-phase, (B) a phase-shifted, and (C) anti-phase relations between the
clusters. The dynamics within the clusters is completely synchronized:
$\varphi_{1}=\varphi_{3}$ and $\varphi_{2}=\varphi_{4}$. These three states
are exemplary, and different phase shifts are obtained from different initial
conditions. In spite of the phase shift, all pendulums are (mean) frequency
synchronized (see the third column of Fig. 4), since they follow the same
motion according to Eq. (20), only phase-shifted. The fourth column of Fig. 4
illustrates the phase-shift between the clusters and their periodic motion
(orange curve). The black point shows the Poincaré map defined by
$\varphi_{1}$ = 0 and $\dot{\varphi}_{1}>0$.

As for the order parameters $r$, Fig. 4(A) presents the trial leading to a
relatively high order parameter close to 1 (almost complete synchronization)
in Fig. 3(B); while Fig. 4(C) represents the trial which falls into the left
side (the inter-group anti-phase synchronization) of the order parameter
distribution in Fig. 3(B). The more trials one draws from random initial
conditions, the more likely one can fill the gap regarding the order parameter
between the inter-group anti-phase synchronization and complete
synchronization to generate EM.

Figure 6:  Dynamical patterns of three types of isolated attractors with
small-amplitude oscillations for 4 coupled pendulums with cross-coupling
structure (13). (A)-(C): SI, IS, and SS synchronization patterns respectively,
(see Table 2 for the pattern explanations). Information along the columns: a:
time series for variables $\theta$, $\varphi_{1}$, $\varphi_{2}$,
$\varphi_{3}$, and $\varphi_{4}$, b: phase-time plots, c: mean frequencies
(not available if a pendulum converges to 0), d: projections on the
$(\varphi_{2},\varphi_{3})$-plane (orange line) and Poincaré maps (black dots,
defined by $\varphi_{2}=0$ and $\dot{\varphi}_{2}>0$, $\varphi_{1}=0$ and
$\dot{\varphi}_{1}>0$, and $\varphi_{1}=0$ and $\dot{\varphi}_{1}>0$, for (A),
(B), and (C), respectively). Parameters are fixed as in Table 1 with
$\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and
$\alpha_{4}=2\pi$ for the 4-coupled clocks ($N=4$).

###  IV.2 Isolated attractors

In addition to infinitely many stable II states from EM described above,
system (13) possesses coexisting isolated attractors corresponding to other
synchronization patterns. These states can also be treated analytically and
numerically in more detail. However, since the main focus of this work is the
EM phenomenon, we consider here exemplary only isolated attractors
corresponding to IA patterns.

####  IV.2.1 Three types of isolated attractors with large-amplitude
oscillations

##### Theoretical analysis of IA patterns.

The IA solutions are characterized by the phase relations
$\varphi_{1}(t)=\varphi_{3}(t)=\psi_{1}$ and
$\varphi_{2}(t)=-\varphi_{4}(t)=\psi_{2}$. Substituting this into system (13),
we obtain:

$\displaystyle\begin{array}[]{ll}\displaystyle(B_{0}+4mr^{2})\ddot{\theta}+k_{\theta}\theta+c_{\theta}\dot{\theta}+mrl[\ddot{\psi}_{2}\sin(\theta-\psi_{2})-\dot{\psi}^{2}_{2}\cos(\theta-\psi_{2})]+mrl[\ddot{\psi}_{2}\sin(\theta+\psi_{2})-\dot{\psi}^{2}_{2}\cos(\theta+\psi_{2})]+\Delta
V_{\theta}=0,\\\
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}+mrl[-\ddot{\theta}\cos(\theta-\psi_{1})+\dot{\theta}^{2}\sin(\theta-\psi_{1})]+\Delta
V_{\varphi_{1}}=M_{E_{1}},\\\
ml^{2}\ddot{\psi}_{2}+mgl\sin{\psi_{2}}+c_{\varphi}{\dot{\psi}_{2}}+mrl[\ddot{\theta}\sin(\theta-\psi_{2})+\dot{\theta}^{2}\cos(\theta-\psi_{2})]+\Delta
V_{\varphi_{2}}=M_{E_{2}},\\\
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}+mrl[\ddot{\theta}\cos(\theta-\psi_{1})-\dot{\theta}^{2}\sin(\theta-\psi_{1})]+\Delta
V_{\varphi_{3}}=M_{E_{3}},\\\
-ml^{2}\ddot{\psi}_{2}-mgl\sin{\psi_{2}}-c_{\varphi}{\dot{\psi}_{2}}+mrl[-\ddot{\theta}\sin(\theta+\psi_{2})-\dot{\theta}^{2}\cos(\theta+\psi_{2})]+\Delta
V_{\varphi_{4}}=M_{E_{4}}.\\\ \end{array}$ (29)

In the following, we introduce the new dimensionless parameter
$\varepsilon=\frac{mrl}{B_{0}+4mr^{2}}$. For the chosen setup as in Table 1,
we have $\varepsilon=0.02726$, i.e., $\varepsilon$ is small, and we will
employ it in our analysis. In fact, the smallness of this parameter is one of
the reasons for the emergence of IA patterns.

Defining further $\bar{k}_{\theta}=\frac{k_{\theta}}{B_{0}+4mr^{2}}=3.7301$,
$\bar{c}_{\theta}=\frac{c_{\theta}}{B_{0}+4mr^{2}}=0.02407$,
$F(\theta,~{}\psi_{2})=-[\ddot{\psi}_{2}\sin(\theta-\psi_{2})-\dot{\psi}^{2}_{2}\cos(\theta-\psi_{2})]+[\ddot{\psi}_{2}\sin(\theta+\psi_{2})-\dot{\psi}^{2}_{2}\cos(\theta+\psi_{2})]-\frac{\Delta
V_{\theta}}{mrl}$ and $\theta=\varepsilon\psi$, the system (29) can be
rewritten in the following form:

$\begin{array}[]{ll}\displaystyle\ddot{\psi}+\bar{c}_{\theta}\dot{\psi}+\bar{k}_{\theta}\psi=F(\varepsilon\psi,~{}\psi_{2}),\\\
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}+\varepsilon
mrl[-\ddot{\psi}\cos(\varepsilon\psi-\psi_{1})+\varepsilon\dot{\psi}^{2}\sin(\varepsilon\psi-\psi_{1})]+\Delta
V_{\varphi_{1}}=M_{E_{1}},\\\
ml^{2}\ddot{\psi}_{2}+mgl\sin{\psi_{2}}+c_{\varphi}{\dot{\psi}_{2}}+\varepsilon
mrl[\ddot{\psi}\sin(\varepsilon\psi-\psi_{2})+\varepsilon\dot{\psi}^{2}\cos(\varepsilon\psi-\psi_{2})]+\Delta
V_{\varphi_{2}}=M_{E_{2}},\\\
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}+\varepsilon
mrl[\ddot{\psi}\cos(\varepsilon{\psi}-\psi_{1})-\varepsilon\dot{\psi}^{2}\sin(\varepsilon{\psi}-\psi_{1})]+\Delta
V_{\varphi_{3}}=M_{E_{3}},\\\
-ml^{2}\ddot{\psi}_{2}-mgl\sin{\psi_{2}}-c_{\varphi}{\dot{\psi}_{2}}+\varepsilon
mrl[-\ddot{\psi}\sin(\varepsilon{\psi}+\psi_{2})-\varepsilon\dot{\psi}^{2}\cos(\varepsilon{\psi}+\psi_{2})]+\Delta
V_{\varphi_{4}}=M_{E_{4}},\end{array}$ (30)

where

$\begin{array}[]{ll}\displaystyle\Delta
V_{\theta}=4lrk_{\varphi}\left(1-\frac{2r}{\sqrt{4r^{2}+2r^{2}(1-\cos(2\psi_{2}))-8lr\sin\psi_{2}\cos(\varepsilon\psi)}}\right)\sin\psi_{2}\sin(\varepsilon\psi),\\\
\Delta V_{\varphi_{1}}=\Delta V_{\varphi_{3}}=0,\\\ \Delta
V_{\varphi_{2}}=k_{\varphi}l\left(1-\frac{2r}{\hat{s}_{24}}\right)[l\sin(2\psi_{2})-2r\cos(\varepsilon\psi-\psi_{2})],\\\
\Delta
V_{\varphi_{4}}=k_{\varphi}l\left(1-\frac{2r}{\hat{s}_{24}}\right)[-l\sin(2\psi_{2})+2r\cos(\varepsilon\psi+\psi_{2})],\end{array}$
(31)

and the values for the distances satisfy:

$\begin{array}[]{ll}\displaystyle s_{13}=s_{31}=s_{24}=s_{42}=2r,\\\
\hat{s}_{13}=\hat{s}_{31}=2r,\\\
\hat{s}_{24}=\hat{s}_{24}=\sqrt{4r^{2}+2r^{2}(1-\cos(2\psi_{2}))-8lr\sin\psi_{2}\cos(\varepsilon\psi)}.\end{array}$
(32)

In the zeroth-order in $\varepsilon$, system (30) is reduced to

$\displaystyle\ddot{\psi}+\bar{c}_{\theta}\dot{\psi}+\bar{k}_{\theta}\psi=F(0,~{}\psi_{2}),$
(33a) $\displaystyle
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}=M_{E_{1}},$
(33b) $\displaystyle
ml^{2}\ddot{\psi}_{2}+mgl\sin{\psi_{2}}+c_{\varphi}{\dot{\psi}_{2}}+\Delta
V_{\varphi_{2}}=M_{E_{2}},$ (33c) $\displaystyle
ml^{2}\ddot{\psi}_{1}+mgl\sin{\psi_{1}}+c_{\varphi}{\dot{\psi}_{1}}=M_{E_{1}},$
(33d) $\displaystyle-
ml^{2}\ddot{\psi}_{2}-mgl\sin{\psi_{2}}-c_{\varphi}{\dot{\psi}_{2}}-\Delta
V_{\varphi_{2}}=-M_{E_{2}},$ (33e)

from which one can see that Eq. (33b) is equivalent to Eq. (33d) and Eq. (33c)
to Eq. (33e). Hence, in this approximation, the subspace of the IA solutions:

$\varphi_{1}(t)=\varphi_{3}(t)=\psi_{1},\quad\varphi_{2}(t)=-\varphi_{4}(t)=\psi_{2},\quad\theta(t)=0.$
(34)

is invariant. For nonzero but small $\varepsilon$, we observe the perturbed
solutions:

$\varphi_{1}(t)\approx\varphi_{3}(t)=\psi_{1},\quad\varphi_{2}(t)\approx-\varphi_{4}(t)=\psi_{2}.$
(35)

Summarizing, the existence of IA patterns can be exactly proven for the limit
$\varepsilon=0$. Since small $\varepsilon$ is a regular perturbation of system
(30), all asymptotically stable periodic attractors in this system will be
only slightly perturbed by small $\varepsilon$-order terms, and one can
observe patterns close to IA.

##### Numerical study of IA, AI, and AA patterns.

The patterns IA, AI, and AA correspond to single lines of the order parameter
distribution in Fig. 3(C)-(E) and, hence, to isolated attractors in the phase
space. Figure 5 reports one example for each of these three patterns. All of
them exhibit partial synchronization with quasiperiodic dynamics.

The phase-time plots in the column b of Fig. 5(A) show that the 1st and 3rd
pendulums are in-phase and the 2nd and 4th are anti-phase with
$\varphi_{1}=\varphi_{3}$ and $\varphi_{2}=-\varphi_{4}$; this is also
confirmed by the analytical solutions (34) and (35). Inside the coupled
groups, the 1st and 3rd or 2nd and 4th pendulums share the same mean
frequency, leading to the multifrequency-clusters Berner  _et al._ (2019);
Omelchenko  _et al._ (2011). Furthermore, the phase space projection on the
plane $(\varphi_{1},~{}\varphi_{2})$ (orange curve) and the corresponding
Poincaré map (black points, when $\varphi_{4}$ = 0 and $\dot{\varphi}_{4}>0$)
indicates that the motion is quasiperiodic. In summary, Fig. 5 reports
different partially synchronous behaviors corresponding to isolated attractors
with large-amplitude osillations: (A) IA, (B) AI, and (C) AA patterns. All
these patterns coexist with EM contributing to a complex multistability
scenario for the cross coupling structure (Fig. 1(G)).

####  IV.2.2 Three types of isolated attractors with small-amplitude
oscillations

The above-mentioned isolated attractors and EM states possess "large-
amplitude" oscillations in the sense that their oscillation amplitude exceeds
the threshold of the escapement mechanism $\varepsilon=5.0^{\circ}$. This
guarantees an inflow of energy into the system and the emergence of stable
self-sustained oscillations of all clocks.

In this section, we investigate the case when some (or all) of the clocks in
system (13) do not reach this threshold. As a result, the SI, IS, and SS
synchronization patterns appear (see Table 2). Figure 6 illustrates the
numerically observed SI, IS, and SS patterns. For example, in Fig. 6(B), both
the time series and phase-time plots exhibit the in-phase synchronization
between $\varphi_{1}$ and $\varphi_{3}$, while $\varphi_{2}$ and $\varphi_{4}$
converge to 0 and stop oscillating. The periodic motions of $\varphi_{1}$ and
$\varphi_{3}$ are illustrated by the phase trajectory (the orange line) and
corresponding Poincaré map (black points, defined by $\varphi_{1}=0$ and
$\dot{\varphi}_{1}>0$) in the column d of Fig. 6(B). Summarizing Fig. 6, it
shows the coexistence of three types of isolated attractors with small-
amplitude oscillations, induced by the discontinuity of the escapement
mechanism.

Figure 7:  One-dimensional dependence of the order parameter on initial
conditions for 4 coupled pendulums with cross-coupling structure (13). The
graph reveals the basins of different attractors along a one-dimensional line
in the phase space. We fix [$\theta^{0}=0.01$,
$\varphi_{1}^{0}=\frac{\pi}{4}$, $\varphi_{3}^{0}=\varphi_{1}^{0}+0.001$,
$\varphi_{4}^{0}=\varphi_{2}^{0}+0.001$,
$\dot{\theta}^{0}=\dot{\varphi}_{1}^{0}=\dot{\varphi}_{2}^{0}=\dot{\varphi}_{3}^{0}=\dot{\varphi}_{4}^{0}=0$],
and vary $\varphi_{2}^{0}$ to obtain 100 discretized values evenly distributed
in the interval $[-\frac{\pi}{4},\frac{\pi}{4}]$. The change of order
parameter with respect to $\varphi_{2}^{0}$ indicates how synchronization
patterns depend on initial conditions, including complete synchronization for
$\varphi_{2}^{0}=\pi/4$. Plateaus correspond to isolated attractors,
continuously changing parts to EM states, and jumps to the basin boundaries.
Parameters are fixed as in Table 1 with $\alpha_{1}=\frac{\pi}{2}$,
$\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and $\alpha_{4}=2\pi$ for the
4-coupled clocks ($N=4$).  Figure 8:  Two-dimensional basins of attraction of
4 coupled pendulums with cross-coupling structure (13). The initial conditions
are chosen as [$\theta^{0}=0.01$,
$\varphi_{1}^{0}=\varphi_{2}^{0}=\frac{\pi}{4}$,
$\varphi_{3}^{0}=\varphi_{1}^{0}-A$, $\varphi_{4}^{0}=\varphi_{2}^{0}-B$,
$\dot{\theta}^{0}=\dot{\varphi}_{1}^{0}=\dot{\varphi}_{2}^{0}=\dot{\varphi}_{3}^{0}=\dot{\varphi}_{4}^{0}=0$],
with $A$ and $B$ changing in the interval $[0,\frac{\pi}{2}]$. The uniform
discretization on the grid $100\times 100$ is used. Four regions in the
bottom-left, top-left, bottom-right, and top-right correspond to the II, IA,
AI, and AA phase patterns, respectively. The bottom-left part with the visible
color gradient correspond to EM. Parameters are fixed as in Table 1 with
$\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and
$\alpha_{4}=2\pi$ for the 4-coupled clocks ($N=4$).  Figure 9:  Two-
dimensional basins of the attractors corresponding to the mixed-mode dynamics
with one cluster staying silent and another oscillating. The model for 4
coupled pendulums with cross-coupling structure (13) are considered. Initial
conditions are chosen as follows: [$\theta^{0}=0.01$,
$\varphi_{1}^{0}=\varphi_{2}^{0}=0$, $\varphi_{3}^{0}=A$, $\varphi_{4}^{0}=B$,
$\dot{\theta}^{0}=\dot{\varphi}_{1}^{0}=\dot{\varphi}_{2}^{0}=\dot{\varphi}_{3}^{0}=\dot{\varphi}_{4}^{0}=0$],
where $A$ and $B$ are taken from the interval $[-\frac{\pi}{4},\frac{\pi}{4}]$
discretized by $100\times 100$ evenly sampled points. The white region
corresponds to the trivial SS pattern (equilibrium at the origin); light green
to SI and IS patterns. Parameters are fixed as in Table 1 with
$\alpha_{1}=\frac{\pi}{2}$, $\alpha_{2}=\pi$, $\alpha_{3}=\frac{3\pi}{2}$ and
$\alpha_{4}=2\pi$ for the 4-coupled clocks ($N=4$).

##  V Basins of attractions

Having clarified the different EM-related asymptotic phase patterns in 4
coupled clocks given by system (13) in Sec. IV, we discuss here their basins
of attractions. Instead of randomly choosing initial conditions for each
pendulum, we fix some of them while the rest are initialized with discretized
values distributed evenly in the given intervals. The corresponding
synchronization states are estimated using order parameter $r$ from Eq. (5),
as it effectively identifies different attractors, including those belonging
to EM. Numerical results are summarized in Figs. 7, 8, and 9.

Figure 7 shows the dependence of the order parameter $r$ on the initial
conditions. For this we fix [$\theta^{0}=0.01$,
$\varphi_{1}^{0}=\frac{\pi}{4}$, $\varphi_{3}^{0}=\varphi_{1}^{0}+0.001$,
$\varphi_{4}^{0}=\varphi_{2}^{0}+0.001$,
$\dot{\theta}^{0}=\dot{\varphi}_{1}^{0}=\dot{\varphi}_{2}^{0}=\dot{\varphi}_{3}^{0}=\dot{\varphi}_{4}^{0}=0$],
and initialize $\varphi_{2}^{0}$ with 100 discretized values evenly
distributed in the interval $[-\frac{\pi}{4},\frac{\pi}{4}]$. We can clearly
observe both the isolated attractors and a part of the EM regime. Isolated
attractors correspond to the flat segments. For example, when
$\varphi_{2}^{0}$ is around 0 in Fig. (7), in spite of different initial
values of $\varphi_{2}^{0}$, trials in this segment part have almost the same
$r$. The continuously changing parts of $r$ in Fig. (7) are in line with the
EM family of II states. The abrupt "jumps" therefore represent boundaries
between either the basins of the isolated attractors or between the isolated
attractors and the EM family.

Figure 8 shows a two-dimensional basin of attraction, where we fix
[$\theta^{0}=0.01$, $\varphi_{1}^{0}=\varphi_{2}^{0}=\frac{\pi}{4}$,
$\varphi_{3}^{0}=\varphi_{1}^{0}-A$, $\varphi_{4}^{0}=\varphi_{2}^{0}-B$,
$\dot{\theta}^{0}=\dot{\varphi}_{1}^{0}=\dot{\varphi}_{2}^{0}=\dot{\varphi}_{3}^{0}=\dot{\varphi}_{4}^{0}=0$],
and vary $A$ and $B$ in the interval $[0,\frac{\pi}{2}]$. The grid of
$100\times 100$ points is used. One can observe four parts corresponding to
II, IA, AI, and AA phase patterns. In particular, the bottom-left part with
the non-constant dependence of the order parameter $r$ in Fig. 8 corresponds
to a subset of an infinite number of stable II states from EM. The other three
parts with constant colors are related to isolated attractors of AI, IA, and
AA patterns. With changing $A$ and $B$, the system moves from complete
synchronization to inter-group anti-phase synchronization (see Fig. 4), and
further to the multifrequency-cluster state (see Fig. 5).

To visualize the basin of the attractors with small-amplitude oscillations, we
analyse a set of initial conditions close to the origin. Figure 9 shows the
corresponding basins of attraction, with the initial conditions
[$\theta^{0}=0.01$, $\varphi_{1}^{0}=\varphi_{2}^{0}=0$, $\varphi_{3}^{0}=A$,
$\varphi_{4}^{0}=B$,
$\dot{\theta}^{0}=\dot{\varphi}_{1}^{0}=\dot{\varphi}_{2}^{0}=\dot{\varphi}_{3}^{0}=\dot{\varphi}_{4}^{0}=0$],
where $A$ and $B$ vary in the interval $[-\frac{\pi}{4},\frac{\pi}{4}]$. In
addition to II patterns from EM set, we obtain here three new regions
corresponding to the patterns SI, IS, and SS, for which some of the clocks are
not oscillating. The basin of the trivial solution SS with all clocks silent
is observed in the central part (white); see Fig. 6(C) for the illustration of
the pattern. The SI and IS patterns are mixed-mode oscillations, with one
cluster silent and another cluster oscillating, corresponding to the light
green basin of attraction.

##  VI Conclusions

In summary, we investigate how different coupling topology affects the
collective dynamics in coupled clocks. The considered model includes global as
well as local couplings represented by the rotating support disc and springs,
respectively. Also, the model contains a discontinuity due to the escapement
clock mechanism. We focus on the model of 4 coupled clocks where surprisingly
an EM is observed between patterns with different synchronization levels. The
EM phenomenon reveals the coexistence of infinitely many stable asymptotic
states. The dependence on initial conditions is clarified using the analysis
of the basins of attractions. The main conclusions based on our work are as
follows:

  * •

For both three coupled and four coupled clocks, we use Monte Carlo sampling
and find that the symmetric coupling structure can increase the dynamical
complexity. This can lead to diverse synchronization patterns and attractors.

  * •

We observe the emergence of EM for the case of four coupled clocks (see Fig.
1(G) and model (13)). This phenomenon is solely induced by the cross-coupling
topological structure, and it is stable against variations of the system
parameters, as far as the clocks remain identical.

  * •

We show both analytically and numerically that the emergence of EM is closely
related to the II synchronization pattern (see Table 2), where the system
splits into two antipodal clusters such that the clocks within the clusters
are fully synchronized, but the intra-cluster dynamics can be shifted by an
arbitrary phase. Moreover, other three types of isolated attractors with
large-amplitude oscillations coexist with the infinite family of states from
the EM family. They correspond to the IA, AI, and AA synchronization patterns.

  * •

We further uncover the effect of the discontinuity of the system induced by
the escapement mechanism. It induces the emergence and coexistence of further
three types of isolated attractors with small-amplitude oscillations. These
states correspond to SI, IS, and SS patterns (see Table 2). In particular, the
IS and SI patterns are mixed states Ebrahimzadeh  _et al._ (2022) where two
clocks are oscillating and the other two stay silent.

As observed in our work, the emergence of EM is caused by a particularly
designed symmetric coupling structure, rather than by introducing additional
quantities into the coupling design. The inclusion of more coupled clocks with
different topological coupling structures is an open way to clarify the
emergence of EM or chimera states of large coupled systems. Another possible
generalization is to include adaptivity in the coupling scheme.

## Acknowledgments

Z.S. and Y.-R.L. was funded by the China Scholarship Council (CSC)
scholarship. J.K. was supported by the Federal Ministry of Education and
Research (BMBF) grant No. 01LP1902J (climXtreme). S.Y. was supported by the
German Research Foundation DFG, Project No. 411803875.

## Author Contributions

Z.S. and Y.S. contributed equally to this work.

## Data Availability

The code used for this work is available online here.

## References

  * Newman (2018) M. Newman,  _Networks_ (Oxford University Press, 2018).
  * Yanchuk  _et al._ (2021) S. Yanchuk, A. C. Roque, E. E. N. Macau, and J. Kurths, The European Physical Journal Special Topics 230, 2711 (2021).
  * Blekhman (1988) I. I. Blekhman,  _Synchronization in Science and Technology_ (ASME Press, New York, NY, USA, 1988).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports 469, 93 (2008).
  * Newman (2003) M. E. J. Newman, SIAM Review 45, 167 (2003), 0303516 .
  * Gross and Blasius (2008) T. Gross and B. Blasius, Journal of The Royal Society Interface 5, 259 (2008), 0709.1858 .
  * Yanchuk  _et al._ (2022) S. Yanchuk, V. Jirsa, J. Cabral, O. Popovych, and A. Torcini,  _From Structure to Function in Neuronal Networks: Effects of Adaptation, Time-Delays, and Noise_, Frontiers Research Topics (Frontiers Media SA, Lausanne, 2022).
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_, Cambridge Nonlinear Science Series (Cambridge University Press, 2001).
  * Osipov  _et al._ (2007) G. V. Osipov, J. Kurths, and C. Zhou,  _Synchronization in Oscillatory Networks_ (Springer Science & Business Media, 2007).
  * Boda  _et al._ (2013) S. Boda, S. Ujvári, A. Tunyagi, and Z. Néda, Eur. J. Phys. 34, 1451 (2013).
  * Chhabria  _et al._ (2018) S. Chhabria, K. A. Blaha, F. Della Rossa, and F. Sorrentino, Chaos 28, 111102 (2018).
  * Ikeguchi and Shimada (2019) T. Ikeguchi and Y. Shimada, in  _Proceedings of the 5th International Conference on Applications in Nonlinear Dynamics_, Understanding Complex Systems, edited by V. In, P. Longhini, and A. Palacios (Springer International Publishing, Cham, 2019) pp. 141–152.
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
  * Dahms  _et al._ (2012) T. Dahms, J. Lehnert, and E. Schöll, Phys. Rev. E 86, 016202 (2012).
  * Lücken and Yanchuk (2012) L. Lücken and S. Yanchuk, Physica D: Nonlinear Phenomena 241, 350 (2012), 1103.6154 .
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004).
  * Omelchenko  _et al._ (2011) I. Omelchenko, Y. Maistrenko, P. Hövel, and E. Schöll, Phys. Rev. Lett. 106, 234102 (2011).
  * Kasatkin  _et al._ (2017) D. V. Kasatkin, S. Yanchuk, E. Schöll, and V. I. Nekorkin, Phys. Rev. E 96, 062211 (2017).
  * Soriano  _et al._ (2013) M. C. Soriano, J. García-Ojalvo, C. R. Mirasso, and I. Fischer, Rev. Mod. Phys. 85, 421 (2013).
  * Tinsley  _et al._ (2012) M. R. Tinsley, S. Nkomo, and K. Showalter, Nature Phys 8, 662 (2012).
  * Hammond  _et al._ (2007) C. Hammond, H. Bergman, and P. Brown, Trends in Neurosciences July INMED/TINS Special Issue—Physiogenic and Pathogenic Oscillations: The Beauty and the Beast, 30, 357 (2007).
  * Blasius  _et al._ (1999) B. Blasius, A. Huppert, and L. Stone, Nature 399, 354 (1999).
  * Tsonis and Roebber (2004) A. A. Tsonis and P. J. Roebber, Physica A: Statistical Mechanics and its Applications 333, 497 (2004).
  * Feudel (2008) U. Feudel, Int. J. Bifurcation Chaos 18, 1607 (2008).
  * Pisarchik and Feudel (2014) A. N. Pisarchik and U. Feudel, Physics Reports Control of Multistability, 540, 167 (2014).
  * Dudkowski  _et al._ (2022) D. Dudkowski, K. Czołczyński, and T. Kapitaniak, Mechanical Systems and Signal Processing 166, 108446 (2022).
  * Arecchi (1991) F. T. Arecchi, Chaos 1, 357 (1991).
  * Foss  _et al._ (1996) J. Foss, A. Longtin, B. Mensour, and J. Milton, Phys. Rev. Lett. 76, 708 (1996).
  * Power and Kleeman (1993) S. B. Power and R. Kleeman, Journal of Physical Oceanography 23, 1670 (1993).
  * Masoller (2002) C. Masoller, Phys. Rev. Lett. 88, 034102 (2002).
  * Borresen and Lynch (2002) J. Borresen and S. Lynch, Int. J. Bifurcation Chaos 12, 129 (2002).
  * Feudel  _et al._ (1996) U. Feudel, C. Grebogi, B. R. Hunt, and J. A. Yorke, Phys. Rev. E 54, 71 (1996).
  * Yanchuk and Wolfrum (2010) S. Yanchuk and M. Wolfrum, SIAM Journal on Applied Dynamical Systems 9, 519 (2010).
  * Balanov  _et al._ (2005) A. G. Balanov, N. B. Janson, and E. Schöll, Phys. Rev. E 71, 016222 (2005).
  * Feudel  _et al._ (1998) U. Feudel, C. Grebogi, L. Poon, and J. A. Yorke, Chaos, Solitons & Fractals 9, 171 (1998).
  * Hens  _et al._ (2015) C. Hens, S. K. Dana, and U. Feudel, Chaos 25, 053112 (2015).
  * Sun  _et al._ (1999) H. Sun, S. K. Scott, and K. Showalter, Phys. Rev. E 60, 3876 (1999).
  * Ngonghala  _et al._ (2011) C. N. Ngonghala, U. Feudel, and K. Showalter, Phys. Rev. E 83, 056206 (2011).
  * Hens  _et al._ (2012) C. R. Hens, R. Banerjee, U. Feudel, and S. K. Dana, Phys. Rev. E 85, 035202 (2012).
  * Pal  _et al._ (2014) S. Pal, B. Sahoo, and S. Poria, Phys. Scr. 89, 045202 (2014).
  * Politi  _et al._ (1986) A. Politi, G. L. Oppo, and R. Badii, Physical Review A 33, 4055 (1986).
  * Pikovsky and Rosenau (2006) A. Pikovsky and P. Rosenau, Physica D: Nonlinear Phenomena 218, 56 (2006).
  * Lamb and Roberts (1998) J. S. Lamb and J. A. Roberts, Physica D: Nonlinear Phenomena 112, 1 (1998).
  * Ashwin  _et al._ (2016) P. Ashwin, C. Bick, and O. Burylko, Frontiers in Applied Mathematics and Statistics 2, 7 (2016).
  * Burylko  _et al._ (2018) O. Burylko, A. Mielke, M. Wolfrum, and S. Yanchuk, SIAM Journal on Applied Dynamical Systems 17, 2076 (2018).
  * Strogatz  _et al._ (2005) S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott, Nature 438, 43 (2005).
  * Dudkowski  _et al._ (2020a) D. Dudkowski, J. Wojewoda, K. Czołczyński, and T. Kapitaniak, Chaos 30, 011102 (2020a).
  * Dudkowski  _et al._ (2020b) D. Dudkowski, P. Jaros, K. Czołczyński, and T. Kapitaniak, Nonlinear Dyn 102, 1541 (2020b).
  * Czołczyński  _et al._ (2011) K. Czołczyński, P. Perlikowski, A. Stefański, and T. Kapitaniak, Chaos 21, 023129 (2011).
  * Kapitaniak  _et al._ (2012) M. Kapitaniak, K. Czolczynski, P. Perlikowski, A. Stefanski, and T. Kapitaniak, Physics Reports Synchronization of Clocks, 517, 1 (2012).
  * Berner  _et al._ (2019) R. Berner, J. Fialkowski, D. V. Kasatkin, V. I. Nekorkin, S. Yanchuk, and E. Schöll, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 103134 (2019).
  * Ebrahimzadeh  _et al._ (2022) P. Ebrahimzadeh, M. Schiek, and Y. Maistrenko, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 103118 (2022), 2206.06830v1 .


# 0.8% Nyquist computational ghost imaging via non-experimental deep learning

Haotian Song  School of Physics, Xi’an Jiaotong University, Xi’an, Shaanxi
710049, China  College of Physics & Astronomy, University of Manchester,
Manchester M13 9PL, UK  Xiaoyu Nie  School of Physics, Xi’an Jiaotong
University, Xi’an, Shaanxi 710049, China  Texas A&M University, College
Station, Texas, 77843, USA  Hairong Su  School of Mathematics and Statistics,
Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China  Hui Chen  Electronic
Materials Research Laboratory, Key Laboratory of the Ministry of Education &
International Center for Dielectric Research, Xi’an Jiaotong University,
Xi’an, 710049, China  Yu Zhou  School of Physics, Xi’an Jiaotong University,
Xi’an, Shaanxi 710049, China  Xingchen Zhao  Texas A&M University, College
Station, Texas, 77843, USA  Tao Peng<EMAIL_ADDRESS>Texas A&M University,
College Station, Texas, 77843, USA  Marlan O. Scully  Texas A&M University,
College Station, Texas, 77843, USA  Baylor University, Waco, 76706, USA
Princeton University, Princeton, NJ 08544, USA

###### Abstract

We present a framework for computational ghost imaging based on deep learning
and customized pink noise speckle patterns. The deep neural network in this
work, which can learn the sensing model and enhance image reconstruction
quality, is trained merely by simulation. To demonstrate the sub-Nyquist level
in our work, the conventional computational ghost imaging results,
reconstructed imaging results using white noise and pink noise via deep
learning are compared under multiple sampling rates at different noise
conditions. We show that the proposed scheme can provide high-quality images
with a sampling rate of 0.8% even when the object is outside the training
dataset, and it is robust to noisy environments. This method is excellent for
various applications, particularly those that require a low sampling rate,
fast reconstruction efficiency, or experience strong noise interference.

††preprint: APS/123-QED

##  I Introduction

Ghost imaging (GI) Pittman  _et al._ (1995); Bennink  _et al._ (2002);
Valencia  _et al._ (2005); Chen  _et al._ (2009) is an innovative method for
measuring the spatial correlations between light beams. With GI, the signal
light field interacts with the object and is collected by a single-pixel
detector, and the reference light field, which does not interact with the
object, falls onto the imaging detector. Therefore, the image information is
not present in either beam alone but only revealed in their correlations.
Computational ghost imaging (CGI) Bromberg  _et al._ (2009); Shapiro (2008)
was proposed to further ameliorate and simplify this framework. In CGI, The
reference arm that records the speckles is replaced by loading pre-generated
patterns directly onto the spatial light modulator or the digital micromirror
device (DMD). The unconventional image is then revealed by correlating the
sequentially recorded intensities at the single-pixel detector with the
corresponding patterns. CGI finds a lot of applications such as wide spectrum
imaging Shrekenhamer  _et al._ (2013); Aspden  _et al._ (2015); Klein  _et
al._ (2019), remote sensing Hardy and Shapiro (2013), and quantum-secured
imaging Clemente  _et al._ (2010).

However, CGI generally requires a large number of samplings to reconstruct a
high-quality image, or the signal would have been submerged under correlation
fluctuations and environmental noise. To suppress the environmental noise and
correlation fluctuations, the required minimum number of sampling is
proportional to the total pixel number of the pattern applied on DMD, i.e.,
the Nyquist sampling limit Cook (1986); Tropp  _et al._ (2009). The image
could have a meager quality with a limited sampling number. This demanding
requirement hindered CGI from fully replacing conventional photography. A
large number of schemes have been proposed to improve CGI’s speed and decrease
the sampling rate (sub-Nyquist). For instance, compressive sensing imaging can
reconstruct images with a relatively low sampling rate by exploiting the
sparsity of the objects Magana-Loaiza  _et al._ (2013); Katz  _et al._ (2009);
Xu  _et al._ (2018); Yi  _et al._ (2019). It nevertheless largely depends on
the sparsity of objects and is sensitive to noise Du  _et al._ (2012).
Orthonormalized noise patterns can be used to suppress the noise and improve
the image’s quality under a limited sampling number Luo  _et al._ (2018); Nie
_et al._ (2020). In particular, the orthonormalized colored noise patterns can
break the Nyquist limit down to $\sim 5\%$ Nie  _et al._ (2020). Fourier and
sequence-ordered Walsh-Hadamard patterns, which are orthogonal to each other
in time or spatial domain, were also applied to the sub-Nyquist imaging Zhang
_et al._ (2015); Wang and Zhao (2016); Zhang  _et al._ (2017). The Russian
doll Sun  _et al._ (2017) and cake-cutting Yu (2019) ordering of Walsh-
Hadamard patterns can minimize the sampling ratio to 5%-10% Nyquist limit.

Recently, the deep learning (DL) technique is employed to identify images He
_et al._ (2020); Li  _et al._ (2020) and improve the quality of images with
the deep neural network (DNN) Lyu  _et al._ (2017); He  _et al._ (2018);
Shimobaba  _et al._ (2018); Wang  _et al._ (2019); Zhai  _et al._ (2019); Wu
_et al._ (2020a); Rizvi  _et al._ (2020); Wu  _et al._ (2020b); Bian  _et al._
(2020). Specifically, computational ghost imaging via deep learning (CGIDL)
has shown a minimum ratio of Nyquist limit down to $\sim 5\%$ He  _et al._
(2018); Wu  _et al._ (2020a). However, such work’s DNNs are trained by
experimental CGI results. Only when the training environment is highly
identical to the environment used for image reconstruction can the DNN be
effective. This limits its universal applications and restricts it to achieve
quick reconstructions. Usually at least thousands of inputs have to be
generated for the training, which would be very time-consuming if conducting
experimental training each time. Some studies have been performed to test the
effectiveness of non-experimental CGI training DNN, the minimum ratios of the
Nyquist limit were up to a few percent Shimobaba  _et al._ (2018); Wang  _et
al._ (2019); Wu  _et al._ (2020b). However, the sampling ratio is much higher
for objects outside of training dataset than those in the training dataset Wu
_et al._ (2020a). Therefore, despite the proliferation of numerous algorithms,
retrieving high-quality images outside of the training group with a meager
Nyquist limit ratio by non-experimental training remains a challenge for the
CGIDL system.

This letter aims to minimize the necessary sampling number further and improve
the imaging quality with the combination of DL and colored noise CGI.
Recently, it has been shown that the synthesized colored noise patterns
possess unique non-zero correlations between neighborhood pixels via amplitude
modulation in the spatial frequency domain Nie  _et al._ (2021); Li  _et al._
(2021). In particular, The pink noise CGI owns positive cross-correlations in
the second-order correlation Nie  _et al._ (2021). It gives a good image
quality under a boisterous environment or pattern distortion when the
traditional CGI method fails. Combining DL with pink noise CGI shows that the
imaging can be retrieved under an extremely low sampling rate ($\sim 0.8\%$).
We also show that we can get training patterns from the simulation without
introducing the environmental noises, i.e., there is no need to get DNN
training with a large number of experimental training inputs. The object used
in the experiment can be independent of the training dataset, which can
largely benefit CGIDL in the real application.

##  II deep learning

###  II.1 Deep neural network

Figure 1: Architecture of DNN. It consists of four convolution layers, one
image input layer, one fully connected layer (yellow), the rectified linear
unit, and the batch normalization layers (red). In the upper line are CGI
results (training inputs) and handwriting ground truths (training labels); In
the bottom line are CGI results from the experiment (test inputs) and CGIDL
results (test outputs) with block style.

Our DNN model, as shown in Fig. 1, uses four convolution layers, one image
input layer, and one fully connected layer. Small $3\times 3$ receptive fields
are applied throughout the whole convolution layers for better performance
Szegedy  _et al._ (2015). Batch normalization layers (BNL), rectified Linear
Unit (ReLU) layers and zero padding are added between each convolution layer.
The BNL is functioned to avoid internal covariate shift during the training
process and speed up the training of DNN Ioffe and Szegedy (2015). The ReLU
layer applies a threshold operation to each element of the inputs Nair and
Hinton (2010). The zero padding part was designed to maintain the
characteristic of input images’ boundaries. To customize the size of training
pictures, both the input and output layers are set to be $54\times 98$. The
solver for training is employed by the Stochastic Gradient Descent with
Momentum Optimizer (SGDMO) to reduce the oscillation via using momentum. The
parameter vector can be updated via equation Eq. (1), which demonstrates the
updating process during the iteration.

$\theta_{\ell+1}=\theta_{\ell}-\alpha\nabla
E\left(\theta_{\ell}\right)+\gamma\left(\theta_{\ell}-\theta_{\ell-1}\right),$
(1)

where $\ell$ is the iteration number, $\alpha$ is the learning rate, $\theta$
is the parameter vector, and $E(\theta)$ is the loss function, mean square
error (MSE). The MSE is defined as

$\mathrm{MSE}=\frac{1}{N_{\mathrm{pixel}}}\sum_{i=1}^{N_{\mathrm{pixel}}}{[\frac{G_{i}-X_{i}}{\langle
G_{(o)}\rangle}]^{2}}.$ (2)

Here, $G$ represents the pixel value of the resulted imaging. $G_{(o)}$
represents pixels that the light ought to be transmitted, i.e., the object
area, while $G_{(b)}$ represents pixels that the light ought to be blocked,
i.e., the background area. $X$ is the ground truth calculated by

${X_{i}}=\begin{cases}\langle G_{(o)}\rangle,&\text{ Transmission = 1}\\\
\langle G_{(b)}\rangle,&\text{ Transmission = 0}\end{cases}$ (3)

The third part on the right hand side of the equation is the feature of SGDMO,
analog to the momentum where $\gamma$ determines the contribution of the
previous gradient step to the current iteration Murphy (2012). Two strategies
are applied to avoid over-fitting of training images. At the end of DNN, a
dropout layer is applied with probability of dropping out input elements being
$0.2$, which is aimed to reduce the connection between convolution layers and
the fully connected layer Srivastava  _et al._ (2014). Meanwhile, we adopted a
step decay schedule for the learning rate. The learning rate dropped from
$10^{-3}$ to $10^{-4}$ after 75 epochs, which constrain the fitting parameters
within a reasonable region. Lower the learning rate could avoid overfitting
significantly with constant maximum epochs.

###  II.2 Network training

Figure 2: The flow chart of CGIDL consists of three parts: (a) training, (b)
test, and (c) experiment. The DNN model is trained with CGI results from
database via simulation. The simulation testing process and experimental
measuments use both the handwriting digits and block style digits. The
experimental part for CGI uses pink noise and white noise speckle patterns,
and their CGI results are ameliorated by trained DNN model.

The proposed CGIDL scheme requires a training process based on pre-prepared
dataset. After training in simulation, it owns ability to reconstruct the
images. We use a set of 10000 digits from the MNIST handwritten digit database
Deng (2012) as training images. All images are resized and normalized to
$54\times 98$ to test the smaller sampling ratio. These training images are
reconstructed by the CGI algorithm. The training images and reconstruction
training images then feed the DNN model as inputs and outputs, respectively,
as shown in Fig. 2(a). The white noise and pink noise speckle patterns are
used separately for the training process, using exactly the same protocol. The
maximum epochs are set as 100, and the training iteration is 31200. The
program is implemented via MATLAB R2019a Update 5 (9.6.0.1174912, 64-bit), and
the DNN is implemented through DL Toolbox. The GPU-chip NVIDIA GTX1050 is used
to accelerate the speed of the computation.

The trained DNN is then tested by simulation and used for retrieving CGI
results in the experiments. In the testing part, the CGI algorithm generates
reconstructed images from testing images with both the MNIST handwritten
digits and block style digits, where the later set is different from images in
the training group. As shown in Fig. 2(b), the trained DNN, fed with
reconstruction testing images, generates CGIDL results. Comparing the
difference between CGIDL and testing images, we could measure the quality of
the trained DNN. Well-performed DNN can be used for retrieving CGI in the
experiment.

The schematic of the experiment is shown in Fig. 2(c). A CW laser is used to
illuminate the DMD, on which the noise patterns are loaded. The pattern
generated by the DMD is then projected onto the object. In our experiment, the
size of the noise patterns is $216\times 392$ DMD pixels ($54\times 98$
independent pixels), in which the independent changeable mirrors count for
$4\times 4$ pixels. Each DMD pixel is $16\mu m\times 16\mu m$ in size.

In the CGI process, the quality of the images is proportional to the sampling
rate, which is the ratio between the number of illumination patterns
$N_{\mathrm{pattern}}$ and $N_{\mathrm{pixel}}$ Ferri  _et al._ (2010); Wang
_et al._ (2015):

$\beta=N_{\mathrm{pattern}}/N_{\mathrm{pixel}}.$ (4)

In the following, We compared the trained network using white noise speckle
patterns (DL white) and pink noise speckle patterns (DL pink), as well as the
conventional CGI (CGI white) in terms of reconstruction performance with
respect to the sampling ratio $\beta$.

##  III Simulation

Figure 3: Simulation results without noise. The upper part used handwriting
digits 1-9 from the training dataset, and the lower part used block style
digits 1-9, which is outside the training dataset. All the simulations are
done at $\beta=5\%$. GT: ground truth.

To test the robustness of our method to different datasets, noise, and its
performs at different sampling rates, we performed a set of simulations. Two
sets of testing images are used in the simulation. One of which is the
handwriting digits 1-9 from the training set, the other is the block style
digits 1-9, which are completely independent of training images. These images
have $28\times 28$ pixels and are resized into $54\times 98$ by widening and
amplification. We started our simulation from the comparison of the CGI white,
DL white and DL pink without noise at $\beta=5\%$, as shown in Fig. 3. The
upper part is with the handwriting digits 1-9, the lower part is with the
block style digits 1-9. Apparently, at this low sampling rate, the traditional
CGI method fails to retrieve the images in both cases. On the other hand, both
DL methods work much better than the traditional CGI. For digits from the
training dataset, both methods work almost equally well. For digits from
outside the training dataset, DL pink works already better than DL white. For
example, the DL white barely can distinguish digits ’3’ and ’8’, but DL pink
can retrieve all the digits images.

In real application, there always exist noise in the measurement. It is
therefore worthwhile checking the performances of different methods under the
influence of noise. We then performed another set of simulations with added
grayscale random noise. The signal-to-noise ratio (SNR) in logarithmic decibel
scale is defined as

$\mathrm{SNR}=10\log\frac{P_{\mathrm{s}}}{P_{\mathrm{b}}},$ (5)

where $P_{\mathrm{s}}$ is the average signal and $P_{\mathrm{b}}$ is the
average noise background. Here the SNR is set to be $4.77\mathrm{dB}$. As
shown in Fig. 4, the upper part is the simulation with digits 2, 3, 5, and 6
from the training dataset, and the lower part is the simulation with digits 2,
3, 5, and 6 from the block style dataset. For both datasets, $\beta$ of 100%,
50%, and 10% are chosen for CGI white, 50%, 5%, and 1% for DL white, 5%, 0.8%,
and 0.5% for DL pink. The image quality is better with the increase of $\beta$
for all cases, as expected. As for the CGI white case, it can only give
marginally visible images when the sampling rate is beyond $50\%$. The DL
white, can retrieve the digits from the training dataset when $\beta=1\%$.
However, for the block style digits, it fails to do so even when $\beta=5\%$.
Unlike the previous case with no noise, there is a significant difference
between objects from the training dataset and outside the training dataset.
Lastly, we note that the DL pink trained network can retrieve the training
dataset when $\beta=0.5\%$. It can also retrieve clear images for the block
style digits at $\beta=0.8\%$. If we compare the black style images at
$\beta=5\%$ for both DL white and DL pink with the no noise case in Fig. 3, it
obvious that the quality of the former is largely affected by the noise, and
the latter is barely affected.

Figure 4: Simulation results of handwriting (top) and block style (bottom)
digits 2, 3, 5, 6 with the SNR of $4.77~{}\mathrm{dB}$. The results of CGI
white are done at $\beta$ of 10%, 50%, and 100%, DL white with $\beta$ of 1%,
5%, and 50%, and DL pink with $\beta$ of 0.5%, 0.8%, and 5%.

##  IV Experiment

Figure 5: Experimental results with the SNR of $14.90~{}\mathrm{dB}$ (upper)
and $4.77~{}\mathrm{dB}$ (lower). Objects are block style digits 2,3, 5, 6.
Different sampling rates are shown for different methods: CGI white are done
at $\beta$ of 10%, 50%, and 100%, DL white with $\beta$ of 1%, 5%, and 50%,
and DL pink with $\beta$ of 0.5%, 0.8%, and 5%.

To further demonstrate the advantage and applicability of CGIDL with pink
noise, we perform experiments with the non-experimental and one-time trained
model. All the experiments are done with digits 2, 3, 5, and 6 with block
style. The block style is chosen to better compare the different behaviors of
all three methods. We manage to start from a relatively low noise level of
$\mathrm{SNR}=14.90~{}\mathrm{dB}$. The results are shown in upper part of Fig
5. We can see at this noise level, the CGI white method barely can distinguish
the images from the noisy background even at $\beta=100\%$. The DL white
trained network, while giving clear images at $\beta=50\%$, fails to fully
image the digits at $\beta=5\%$. This is mainly due to the objects are outside
the training set, reveal one of the shortcomings of the standard DL network.
On the other hand, our DL pink trained network can still give clear results
even when $\beta=0.5\%$.

We then increase the noise level to $\mathrm{SNR}=4.77~{}\mathrm{dB}$, which
is the same as the simulation case so we can have a fair comparison. The
experimental results are shown in the lower part of Fig 5. The CGI white
completely fail to image the digits even at $\beta=100\%$. The DL white
trained network is also largely affected by the noisy environment, and not
able to fully retrieve the images at $\beta=50\%$. On the other hand, the DL
pink method can still image all digits at the sampling rate of $0.8\%$. If we
compare these results to the corresponding low noise case, we can see that the
image qualities do not change much, indicating our trained network is robust
to noise. Also, the results with $\beta=0.8\%$ is better than the standard DL
white network at $\beta=50\%$, which is about two orders higher.

To quantitatively justify the quality of reconstructed block style images, we
compare three evaluating indicators of image quality, i.e., the peak signal to
noise ratio (PSNR), the visibility (VIS), and the correlation coefficient
(CC):

$\displaystyle\mathrm{PSNR}$
$\displaystyle=10\times{\log_{10}[\frac{(2^{k}-1)^{2}}{\mathrm{MSE}}]},$ (6)
$\displaystyle\mathrm{VIS}$ $\displaystyle=\frac{\langle
G_{(\mathrm{o})}\rangle-\langle G_{(\mathrm{b})}\rangle}{\langle
G_{(\mathrm{o})}\rangle+\langle G_{(\mathrm{b})}\rangle},$ (7)
$\displaystyle\mathrm{CC}$
$\displaystyle=\frac{\mathrm{Cov}(G,X)}{\sqrt{\mathrm{Var}(G)\mathrm{Var}(X)}}.$
(8)

Here $\mathrm{MSE}$ is defined in Eq. 2, $\mathrm{Var(\,)}$ is the variance of
its arguments, $\mathrm{Cov(\,)}$ is the covariance of its arguments, $k$ is
the gray level of the image, and in our experiment $k\equiv 8$.

Figure 6: PSNR, VIS, and CC for simulation and experiments of the block style
digits 2, 3, 5, 6 with three categories: CGI white ($\beta$ of 10%, 50%, and
100%), DL white ($\beta$ of 1%, 5%, and 50%), and DL pink ($\beta$ of 0.5%,
0.8%, and 5%).

The results for all cases including simulation without and with noise,
experiment with high and low SNR, are shown in Fig. 6. The $\mathrm{PSNR}$,
$\mathrm{VIS}$, and $\mathrm{CC}$ all indicate that the CGIDL methods are much
better than the traditional CGI method. Indeed, as shown in the simulation
case, the image quality of CGIDL at 5% is already better than CGI at full
sampling rate for all the situations. When we compare the two DL methods, we
see that in general DL pink is much better than DL white, as also suggested
from Figs. 3, 4, and 5. Since the network is trained using $\mathrm{MSE}$ as
the loss function, the $\mathrm{PSNR}$ of simulation without noise at 5% is
very similar for both cases. However, when the noise increases, the
$\mathrm{PSNR}$ of DL white starts to decrease, while the $\mathrm{PSNR}$ of
DL pink does not change much. The $\mathrm{VIS}$ and $\mathrm{CC}$ also have
similar behavior as $\mathrm{PSNR}$. We note here that all three indicators
suggest DL pink works better than the other two methods, in the experimental
results with low SNR, DL pink of 5% sampling rate is already better than the
DL white with 50% sampling rate.

##  V Conclusion

In conclusion, we have demonstrated a deep-learning imaging method with pink
noise patterns. The DNN is trained using only simulation data from the
handwriting dataset. The trained network can then be applied to various
conditions, including objects outside the training set and experiments with
strong noise. We have demonstrated imaging results with extremely low sampling
rate both in simulation and experiments. We have also evaluated the quality of
the images outside the training dataset for both simulation and experimental
results, in terms of $\mathrm{PSNR}$, $\mathrm{VIS}$, and $\mathrm{CC}$.

All results suggest that the DL pink scheme has a great advantage, especially
in the low sampling region. This one-time, noise-robust, and non-experimental
training CGIDL is eligible to be implemented in various situations and has a
wide range of application prospects. The pink noise speckle patterns, trained
DNN with various sampling rates, and their raw encoding programs are
encapsulated and uploaded online 111https://github.com/XJTU-TAMU-CGI/CGIDL.
People who need a quick sampling function on CGIDL can utilize this universal
system to get ameliorated results in other CGIDL systems. Further works can
reach to other imaging and spectroscopy systems by loss function adjustment
and speckle pattern optimization, in order to get spatial, frequency, or time-
resolution. In addition to results amelioration, DL may also have great
potential to generate optimized speckle patterns for a variety of tasks.

## Funding

Air Force Office of Scientific Research (Award No. FA9550-20-1-0366 DEF),
Office of Naval Research (Award No. N00014-20-1-2184), Robert A. Welch
Foundation (Grant No. A-1261), National Science Foundation (Grant No.
PHY-2013771).

## Data availability.

The experimental and simulation data are available upon reasonable request.
The trained DL networks and raw DL training and test codes are uploaded on the
website: https://github.com/XJTU-TAMU-CGI/CGIDL.

## Disclosures

The authors declare no conflicts of interest.

## References

  * Pittman  _et al._ (1995) T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, Optical imaging by means of two-photon quantum entanglement, Physical Review A 52, R3429 (1995).
  * Bennink  _et al._ (2002) R. S. Bennink, S. J. Bentley, and R. W. Boyd, “two-photon” coincidence imaging with a classical source, Physical Review Letters 89, 113601 (2002).
  * Valencia  _et al._ (2005) A. Valencia, G. Scarcelli, M. D’Angelo, and Y. Shih, Two-photon imaging with thermal light, Physical Review Letters 94, 063601 (2005).
  * Chen  _et al._ (2009) X.-H. Chen, Q. Liu, K.-H. Luo, and L.-A. Wu, Lensless ghost imaging with true thermal light, Optics Letters 34, 695 (2009).
  * Bromberg  _et al._ (2009) Y. Bromberg, O. Katz, and Y. Silberberg, Ghost imaging with a single detector, Physical Review A 79, 053840 (2009).
  * Shapiro (2008) J. H. Shapiro, Computational ghost imaging, Physical Review A 78, 061802 (2008).
  * Shrekenhamer  _et al._ (2013) D. Shrekenhamer, C. M. Watts, and W. J. Padilla, Terahertz single pixel imaging with an optically controlled dynamic spatial light modulator, Optics Express 21, 12507 (2013).
  * Aspden  _et al._ (2015) R. S. Aspden, N. R. Gemmell, P. A. Morris, D. S. Tasca, L. Mertens, M. G. Tanner, R. A. Kirkwood, A. Ruggeri, A. Tosi, R. W. Boyd, _et al._ , Photon-sparse microscopy: visible light imaging using infrared illumination, Optica 2, 1049 (2015).
  * Klein  _et al._ (2019) Y. Klein, A. Schori, I. Dolbnya, K. Sawhney, and S. Shwartz, X-ray computational ghost imaging with single-pixel detector, Optics Express 27, 3284 (2019).
  * Hardy and Shapiro (2013) N. D. Hardy and J. H. Shapiro, Computational ghost imaging versus imaging laser radar for three-dimensional imaging, Physical Review A 87, 023820 (2013).
  * Clemente  _et al._ (2010) P. Clemente, V. Durán, E. Tajahuerce, J. Lancis, _et al._ , Optical encryption based on computational ghost imaging, Optics Letters 35, 2391 (2010).
  * Cook (1986) R. L. Cook, Stochastic sampling in computer graphics, ACM Transactions on Graphics (TOG) 5, 51 (1986).
  * Tropp  _et al._ (2009) J. A. Tropp, J. N. Laska, M. F. Duarte, J. K. Romberg, and R. G. Baraniuk, Beyond nyquist: Efficient sampling of sparse bandlimited signals, IEEE transactions on information theory 56, 520 (2009).
  * Magana-Loaiza  _et al._ (2013) O. S. Magana-Loaiza, G. A. Howland, M. Malik, and J. C. Howell, Compressive object tracking using entangled photons, Applied Physics Letters 102 (2013).
  * Katz  _et al._ (2009) O. Katz, Y. Bromberg, and Y. Silberberg, Compressive ghost imaging, Applied Physics Letters 95, 739 (2009).
  * Xu  _et al._ (2018) Z.-H. Xu, W. Chen, J. Penuelas, M. Padgett, and M.-J. Sun, 1000 fps computational ghost imaging using led-based structured illumination, Opt. Express 26, 2427 (2018).
  * Yi  _et al._ (2019) C. Yi, C. Zhengdong, F. Xiang, C. Yubao, and L. Zhenyu, Compressive sensing ghost imaging based on image gradient, Optik 182, 1021 (2019).
  * Du  _et al._ (2012) J. Du, W. Gong, and S. Han, The influence of sparsity property of images on ghost imaging with thermal light, Optics Letters 37, 1067 (2012).
  * Luo  _et al._ (2018) B. Luo, P. Yin, L. Yin, G. Wu, and H. Guo, Orthonormalization method in ghost imaging, Optics express 26, 23093 (2018).
  * Nie  _et al._ (2020) X. Nie, X. Zhao, T. Peng, and M. O. Scully, Sub-nyquist computational ghost imaging with orthonormalized colored noise pattern, arXiv preprint arXiv:2012.07250 (2020).
  * Zhang  _et al._ (2015) Z. Zhang, X. Ma, and J. Zhong, Single-pixel imaging by means of fourier spectrum acquisition, Nature Communications 6, 1 (2015).
  * Wang and Zhao (2016) L. Wang and S. Zhao, Fast reconstructed and high-quality ghost imaging with fast walsh–hadamard transform, Photonics Research 4, 240 (2016).
  * Zhang  _et al._ (2017) Z. Zhang, X. Wang, G. Zheng, and J. Zhong, Hadamard single-pixel imaging versus fourier single-pixel imaging, Opt. Express 25, 19619 (2017).
  * Sun  _et al._ (2017) M.-J. Sun, L.-T. Meng, M. P. Edgar, M. J. Padgett, and N. Radwell, A russian dolls ordering of the hadamard basis for compressive single-pixel imaging, Scientific Reports 7, 1 (2017).
  * Yu (2019) W.-K. Yu, Super sub-nyquist single-pixel imaging by means of cake-cutting hadamard basis sort, Sensors 19, 4122 (2019).
  * He  _et al._ (2020) X. He, S.-M. Zhao, and L. Wang, Handwritten digit recognition based on ghost imaging with deep learning, Chinese Physics B (2020).
  * Li  _et al._ (2020) J. Li, M. Le, J. Wang, W. Zhang, B. Li, and J. Peng, Object identification in computational ghost imaging based on deep learning, Applied Physics B 126, 1 (2020).
  * Lyu  _et al._ (2017) M. Lyu, W. Wang, H. Wang, H. Wang, G. Li, N. Chen, and G. Situ, Deep-learning-based ghost imaging, Scientific Reports 7, 17865 (2017).
  * He  _et al._ (2018) Y. He, G. Wang, G. Dong, S. Zhu, H. Chen, A. Zhang, and Z. Xu, Ghost imaging based on deep learning, Scientific Reports 8, 6469 (2018).
  * Shimobaba  _et al._ (2018) T. Shimobaba, Y. Endo, T. Nishitsuji, T. Takahashi, Y. Nagahama, S. Hasegawa, M. Sano, R. Hirayama, T. Kakue, A. Shiraki, _et al._ , Computational ghost imaging using deep learning, Optics Communications 413, 147 (2018).
  * Wang  _et al._ (2019) F. Wang, H. Wang, H. Wang, G. Li, and G. Situ, Learning from simulation: An end-to-end deep-learning approach for computational ghost imaging, Optics Express 27, 25560 (2019).
  * Zhai  _et al._ (2019) X. Zhai, Z.-d. Cheng, Y. Chen, Z.-y. Liang, Y. Wei, _et al._ , Foveated ghost imaging based on deep learning, Optics Communications 448, 69 (2019).
  * Wu  _et al._ (2020a) H. Wu, R. Wang, G. Zhao, H. Xiao, D. Wang, J. Liang, X. Tian, L. Cheng, and X. Zhang, Sub-nyquist computational ghost imaging with deep learning, Optics Express 28, 3846 (2020a).
  * Rizvi  _et al._ (2020) S. Rizvi, J. Cao, K. Zhang, and Q. Hao, Deepghost: real-time computational ghost imaging via deep learning, Scientific Reports 10, 1 (2020).
  * Wu  _et al._ (2020b) H. Wu, R. Wang, G. Zhao, H. Xiao, J. Liang, D. Wang, X. Tian, L. Cheng, and X. Zhang, Deep-learning denoising computational ghost imaging, Optics and Lasers in Engineering 134, 106183 (2020b).
  * Bian  _et al._ (2020) T. Bian, Y. Yi, J. Hu, Y. Zhang, Y. Wang, and L. Gao, A residual-based deep learning approach for ghost imaging, Scientific Reports 10, 1 (2020).
  * Nie  _et al._ (2021) X. Nie, F. Yang, X. Liu, X. Zhao, R. Nessler, T. Peng, M. S. Zubairy, and M. O. Scully, Noise-robust computational ghost imaging with pink noise speckle patterns, Phys. Rev. A 104, 013513 (2021).
  * Li  _et al._ (2021) Z. Li, X. Nie, F. Yang, X. Liu, D. Liu, X. Dong, X. Zhao, T. Peng, M. S. Zubairy, and M. O. Scully, Sub-rayleigh second-order correlation imaging using spatially distributive colored noise speckle patterns, Optics Express 29, 19621 (2021).
  * Szegedy  _et al._ (2015) C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, Going deeper with convolutions, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2015) pp. 1–9.
  * Ioffe and Szegedy (2015) S. Ioffe and C. Szegedy, Batch normalization: Accelerating deep network training by reducing internal covariate shift, arXiv preprint arXiv:1502.03167 (2015).
  * Nair and Hinton (2010) V. Nair and G. E. Hinton, Rectified linear units improve restricted boltzmann machines, in  _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , ICML’10 (Omnipress, Madison, WI, USA, 2010) pp. 807–814.
  * Murphy (2012) K. P. Murphy,  _Machine learning: a probabilistic perspective_ (MIT press, 2012).
  * Srivastava  _et al._ (2014) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, Dropout: A simple way to prevent neural networks from overfitting, Journal of Machine Learning Research 15, 1929 (2014).
  * Deng (2012) L. Deng, The mnist database of handwritten digit images for machine learning research [best of the web], IEEE Signal Processing Magazine 29, 141 (2012).
  * Ferri  _et al._ (2010) F. Ferri, D. Magatti, L. Lugiato, and A. Gatti, Differential ghost imaging, Physical review letters 104, 253603 (2010).
  * Wang  _et al._ (2015) W. Wang, X. Hu, J. Liu, S. Zhang, J. Suo, and G. Situ, Gerchberg-saxton-like ghost imaging, Optics Express 23, 28416 (2015).
  * Note (1) Https://github.com/XJTU-TAMU-CGI/CGIDL.