ABSTRACT
Title of dissertation: EXPERIMENTAL STUDY OF
ENTROPY GENERATION AND
SYNCHRONIZATION IN NETWORKS
OF COUPLED DYNAMICAL SYSTEMS
Aaron Morgan Hagerstrom
Doctor of Philosophy, 2015
Dissertation directed by: Professor Rajarshi Roy
Department of Physics
This thesis describes work in two areas: unpredictability in a system with
both single photon detection and chaos, and synchronization in networks of coupled
systems.
The unpredictability of physical systems can depend on the scale at which
they are observed. For example, single photons incident on a detector arrive at ran-
dom times, but slow intensity variations can be observed by counting many photons
over large time windows. We describe an experiment in which a weak optical sig-
nal is modulated by feedback from a single-photon detector. We demonstrate that
at low photon rates, the photon arrivals are described by Poisson noise, while at
high photon rates, the intensity of the light shows a deterministic chaotic modula-
tion. Furthermore, we show that measurements of the entropy rate can resolve both
noise and chaos, and describe the relevance of this observation to random number
generation.
We also describe an experimental system that allows for the study of large
networks of coupled dynamical systems. This system uses a spatial light modulator
(SLM) and a camera in a feedback configuration, and has an optical nonlinearity
due to the relationship between a spatially-dependent phase shift applied by the
modulator, and the intensity measured by the camera. This system is a powerful
platform for studying synchronization in large networks of coupled systems, because
it can iterate many (∼1000) maps in parallel, and can be configured with any desired
coupling matrix. We use this system to study cluster synchronization. It is often
observed that networks synchronize in clusters. We show that the clusters that form
can be predicted based on the symmetries of the network, and their stability can be
analyzed. This analysis is supported by experimental measurements.
The SLM feedback system can also exhibit chimera states. These states were
first seen in models of large populations of coupled phase oscillators. Their defining
feature is the coexistence of large populations of synchronized and unsynchronized
oscillators in spite of symmetrical coupling configurations and homogeneous oscilla-
tors. Our SLM feedback system provided one of the first two experimental examples
of chimera states. We observed chimeras in coupled map lattices. These states evolve
in discrete time, and in this context, “coherence” and “incoherence” refer to smooth
and discontinuous regions of a spatial pattern.
EXPERIMENTAL STUDY OF ENTROPY GENERATION AND
SYNCHRONIZATION IN NETWORKS OF COUPLED
DYNAMICAL SYSTEMS
by
Aaron Morgan Hagerstrom
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2015
Advisory Committee:
Professor Rajarshi Roy, Chair/Advisor
Professor Thomas E. Murphy, Co-Advisor
Dr. Louis M. Pecora
Professor Christopher Jarzynski
Professor Daniel Lathrop
Professor Edward Ott
c© Copyright by
Aaron Morgan Hagerstrom
2015
Table of Contents
List of Figures iv
1 Introduction 1
1.1 Entropy generation in a photon-counting feedback loop . . . . . . . . 2
1.2 Synchronization in networks of coupled dynamical systems . . . . . . 4
1.3 Organization of this thesis . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Photon Counting Optoelectronic Feedback Loop 8
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Motivation and background . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.1 Discrete and analog dynamics . . . . . . . . . . . . . . . . . . 11
2.2.2 Components of the feedback loop . . . . . . . . . . . . . . . . 14
2.3 Digital Implementation on FPGA . . . . . . . . . . . . . . . . . . . . 16
2.3.1 Photon Detection and Time delay . . . . . . . . . . . . . . . . 17
2.3.2 Filter equations . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.1 Time series, probability distributions, and autocorrelation func-
tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4.2 Attractor visualizations . . . . . . . . . . . . . . . . . . . . . . 23
2.4.3 Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 Entropy generation in a photon-driven feedback loop 27
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.1 NIST Recommendations for Entropy Sources . . . . . . . . . . 31
3.2.2 Entropy Generation in Physical Systems . . . . . . . . . . . . 32
3.3 Entropy and Entropy Rate . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 Entropy Calculations . . . . . . . . . . . . . . . . . . . . . . . 38
3.4 Entropy rate of a constant rate Poisson process . . . . . . . . . . . . 40
3.4.1 Entropy rate from inter-arrival intervals . . . . . . . . . . . . 42
3.4.2 Can intensity modulation increase entropy rates? . . . . . . . 43
3.4.3 Counting windows and ε . . . . . . . . . . . . . . . . . . . . . 44
ii
3.5 Entropy rate of chaotic systems . . . . . . . . . . . . . . . . . . . . . 45
3.6 Entropy generation in a photon-driven feedback loop . . . . . . . . . 49
3.7 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.7.1 Broader context . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4 Cluster Synchronization and Symmetries in Networks 55
4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2.1 Coupling – adjacency matrix and Laplacian matrix . . . . . . 57
4.2.2 General model of networks of coupled identical systems . . . . 58
4.3 Symmetries and cluster synchronization . . . . . . . . . . . . . . . . . 60
4.4 Master stability function formalism . . . . . . . . . . . . . . . . . . . 64
4.5 Variational equations in IRR Coordinates . . . . . . . . . . . . . . . . 66
4.6 Calculating the IRR transformation . . . . . . . . . . . . . . . . . . . 69
4.7 Experimental system . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.8 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.8.1 Subgroup decomposition and isolated desynchronization . . . 77
4.9 Generalizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.9.1 Subgroup clusters . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.9.2 Laplacian coupling and merged clusters . . . . . . . . . . . . . 81
4.9.3 Experimental observation . . . . . . . . . . . . . . . . . . . . 83
5 Chimeras in Coupled-Map Lattices 85
5.1 Introduction and outline . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.1.1 Organization of this chapter . . . . . . . . . . . . . . . . . . . 88
5.2 Coupled-map lattices in experiment . . . . . . . . . . . . . . . . . . . 89
5.2.1 Chimeras in 1D and 2D . . . . . . . . . . . . . . . . . . . . . 89
5.2.2 Periodicity quantification . . . . . . . . . . . . . . . . . . . . . 95
5.3 Critical coupling strength . . . . . . . . . . . . . . . . . . . . . . . . 96
5.4 Scaling of the incoherent regions in chimera states with N . . . . . . . 101
5.5 Spatial rescaling and coherent patterns with K=1,2,3... . . . . . . . . 104
5.6 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6 Unanswered questions 107
6.1 Entropy and the photon-counting feedback loop . . . . . . . . . . . . 107
6.2 Synchronization in networks of coupled systems . . . . . . . . . . . . 109
Bibliography 112
iii
List of Figures
2.1 Analog and photon-counting optoelectronic feedback loops. . . . . . . 12
2.2 Time series, probability distributions, and autocorrelation functions . 21
2.3 Poincare´ sections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4 Experimental dependence of growth rate of variance on counting time
window, w. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1 Illustration of discretization in time (τ) and value (ε). . . . . . . . . . 33
3.2 Pattern entropy and entropy rate . . . . . . . . . . . . . . . . . . . . 37
3.3 Entropy of a constant-rate Poisson process with a counting window . 46
3.4 Entropy calculation from a deterministic model of the photon-counting
feedback loop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.5 Entropy rate of experimental time series. . . . . . . . . . . . . . . . . 50
4.1 Symmetries and clusters in a simple network. . . . . . . . . . . . . . 61
4.2 Three randomly generated networks with varying amounts of sym-
metry and associated coupling matrices. . . . . . . . . . . . . . . . . 63
4.3 Experimental configuration. . . . . . . . . . . . . . . . . . . . . . . . 73
4.4 Experimental observation of isolated and intertwined desynchroniza-
tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.5 All of the synchronization patterns of a small network. . . . . . . . . 82
4.6 Experimental demonstration of subgroup and Laplacian clusters . . . 84
5.1 Experimental device and coupling configuration . . . . . . . . . . . . 90
5.2 Parameter space of the 1D and 2D CML systems. . . . . . . . . . . . 92
5.3 Critical coupling and scaling of coherent profiles. . . . . . . . . . . . . 97
5.4 Local map and its twice iterate . . . . . . . . . . . . . . . . . . . . . 99
5.5 Simulations to illustrate scaling with system size of the incoherent
regions for the 1D system. . . . . . . . . . . . . . . . . . . . . . . . . 102
5.6 Coherent profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
iv
Chapter 1: Introduction
Like many important concepts in science, chaos was discovered several times.
Poincare´ encountered chaos while studying the three-body problem in orbital dy-
namics [1] in the early 20th century. Chaos was also experimentally observed in
the sinusoidally-driven van der Pol oscillator 1927, and described mathematically in
this context the 1950s [2, 3]. One of the first scientists to realize the far-reaching
importance of chaos was Lorenz, who discovered chaos while working on a simplified
model of convection in the atmosphere. His paper from 1963, Deterministic Nonpe-
riodic Flow [4] described many ideas that are central to the modern understanding of
chaos, including dissipative dynamics in phase space, the fractal structure of attrac-
tors (which he described in a qualitative way), and sensitivity to initial conditions.
This sensitivity to initial conditions is the property of chaos that allows for random
number generation.
The idea that chaos can be used to generate random numbers is old – older
than the use of the word “chaos” to describe unpredictable nonperiodic motion in
deterministic systems. The word “chaos” has been in use since 1975 [5], although
its meaning has evolved somewhat over the years [6].
1
In 1947 [7, 8] Ulam and von Neumann showed that the logistic map,
xt+1 = 4xt(1− xt), (1.1)
generates one bit of information with every iteration. They introduced the change
of variables xt = sin2(pizt), and showed that Equation 1.1 has an explicit solution.
zt = 2
tz0 mod 1 (1.2)
This solution gives a concrete example of how a deterministic system can generate
information. If the initial condition, z0, is written in binary, then, according to
Equation 1.2, the decimal place is moved to the right one place for each iteration.
The trailing, insignificant digits of z become more significant as time passes. If,
at each iteration, one reads out the binary digit to the right of the decimal place,
one will find a string of random bits that is as unbiased as a fair coin. There are
specific values of z0 for which this is not the case, for example 0.11111111111111...,
and 0.10101010101010...., but almost all real numbers on the interval [0,1] have as
many 0’s as 1’s, and in general, an equal number of occurrences of the bit patterns
of length d. In other words, almost all real numbers have one bit of entropy per
binary digit [9]. By amplifying differences in initial conditions, chaos makes these
insignificant trailing bits of information large enough to measure [10].
1.1 Entropy generation in a photon-counting feedback loop
Since in 2008, chaotic semiconductor lasers have been used to generate random
numbers [11]. These systems are appealing for (at least) two reasons. First, they
2
are fast (on the order of 109 to 1012 bits per second) [12]. Second, random bits
generated by physical processes offer additional cryptographic security compared
to pseudorandom algorithms because their unpredictability can be attributed to a
fundamentally unpredictable physical process.
Not all physical random number generators are based on chaos. Another, much
longer used, physical source of unpredictability is the detection of single photons
[13, 14]. Several commercial devices, which are marketed as “quantum random
number generators” are currently being sold [15], which are based on weak light
sources and photon-counting detectors. One scheme employs a beam splitter and
two detectors, so that whenever a photon arrives, one detector or the other detects
the photon, yielding either a 0 or a 1. Another scheme involves harvesting entropy
from the random time intervals between photon arrivals [16].
In this thesis, we describe an experimental system with unpredictability due
to both chaotic dynamics and discrete photon arrivals. It similar to earlier opto-
electronic feedback loops [17,18]. These feedback loops are known to have complex
dynamics that can be tuned anywhere between fixed points and high-dimensional
chaos, and whose timescales can be tuned over a wide range [19,20]. By incorporat-
ing a photon-counting detector, we are able to observe a transition from shot noise
to chaos with increasing photon rate. Furthermore, we find that the type of unpre-
dictability we observe – shot noise or chaos – depends on our scales of observation
in both time and value.
These results have immediate relevance to physical random number genera-
tion. Usually, physical random number generators are evaluated based on their
3
ability to produce bit streams free of bias and correlation [21,22]. Physical systems
usually have both bias (for example, non-identical detectors) and correlation (dy-
namical correlation times, finite detector bandwidth). In practice, post-processing
is very often employed to improve the statistical properties of a noise source. It
can be difficult to assess whether the entropy in the resulting bit stream should
be attributed to the physical noise source or the post-processing procedure. Fur-
thermore, the designer of a physical random number generator must decide how
frequently to sample their system and to what resolution they should measure it.
So, they should be aware of the effect of coarse graining on the entropy rate of their
system. We show that the scaling of the entropy rate with sampling resolution can
be attributed to either chaos or shot noise depending on the photon count rate and
sampling resolution.
1.2 Synchronization in networks of coupled dynamical systems
Often, networks of coupled dynamical systems will synchronize in clusters. Up
until a few years ago, only anecdotal examples could be analyzed in detail [23–26].
We show here that the symmetries of a network can be used to predict the clusters
that form in an arbitrary network. The master stability function formalism has
proven to be a very useful paradigm for the study of global synchronization in
networks of coupled systems [27]. This technique examines the linear stability of
the globally-synchronous solution. We present here the linear stability analysis of
cluster synchronization, and show that the symmetries of a network can be used to
4
determine a coordinate transformation that vastly simplifies the this analysis. Now,
due to work presented here, it is possible to predict and analyze the stability of a
great number of cluster-synchronous states in a given arbitrary network. We also
present an experimental verification of the predictions made by this technique. In
the experiment, a novel bifurcation called “isolated synchronization” is observed,
where a single cluster looses synchrony while others remain synchronized. This
bifurcation can also be explained in terms of the symmetries of the network.
In larger networks, one is often interested in dynamical behaviors that persist
in the the “thermodynamic limit”, where the number of oscillators becomes infinite.
The Kuramoto model of coupled phase oscillators [28] is particularly appealing be-
cause it can be shown analytically that in the thermodynamic limit a population of
phase oscillators undergoes a phase transition from unsynchronized to synchronized
behavior with increasing coupling strength.
In the early 2000’s, it was realized that in certain coupling configurations,
Kuramoto oscillators could form “chimera states” in which there were coexisting
synchronized and unsynchronized populations [29–31]. It was not until 2012 that
such states were observed in experiment [32,33]. Now, there are several quite diverse
experimental examples of chimera states [32–38]. Here, we describe one of the first
two experimental observations of chimera states.
5
1.3 Organization of this thesis
Chapter 2 describes the photon-counting feedback loop. The operation of
the components in the loop will be described, as well as some details of the digital
implementation of the signal processing and data acquisition. Finally, the transition
between shot noise and deterministic chaos will be described, and illustrated by
examining time series plots, Poincare´ sections, autocorrelation functions, and the
growth of the variance in the number of photon counts with the counting time
window.
Chapter 3 describes the quantification of entropy production in the photon-
counting feedback loop. This chapter also contains brief descriptions of some basic
concepts from information theory, and a couple of different analyses of the entropy
rate of a Poisson process. This material, while not new, provides an illustration
of important concepts and a point of reference for understanding the experimental
results. The main goal of this chapter is to show that the entropy rate of the photon
counting feedback loop reflects chaos at large scales and high photon rates, and shot
noise at small scales and low photon rates.
Chapter 4 describes cluster synchronization. A lot of space will be devoted
to describing the linear stability analysis of the cluster-synchronous state, and how
the symmetries of the network are crucial in predicting with clusters form and
analyzing their stability. We experimentally demonstrate cluster synchrony in a few
networks. In one of these networks, the symmetries are far from obvious, and we
show that computational group theory can be used to analyze cluster synchrony in
6
this network. In another network, the symmetries are obvious by inspection, and
we use this network as an example to show both the diversity of possible patterns
and the method for predicting them and analyzing their stability.
Chapter 5 describes an experimental observation of chimera states. These
chimera states exist in a coupled map lattice, and have some qualitative similarities
and differences compared with chimera states in populations of phase oscillators. A
one-dimensional and a two-dimensional version of the same system are be compared,
and a few basic theoretical results are be given.
Chapter 6 presents a list of questions. These are ideas that arose during our
investigation that I did not have time to pursue.
7
Chapter 2: Photon Counting Optoelectronic Feedback Loop
This chapter is based on work from the following publication
Aaron Morgan Hagerstrom, Thomas Edward Murphy, and Rajarshi Roy. Harvest-
ing entropy and quantifying the transition from noise to chaos in a photon-counting
feedback loop. Proceedings of the National Academy of Sciences 112(30):9258–9263,
2015.
2.1 Overview
This chapter describes an optoelectronic feedback loop in which noise due
to discrete photon arrivals plays an important role in the dynamics. The chapter
proceeds as follows: Section 2.2 describes the motivation for this experiment, and
highlights some previous work on similar experimental systems in order to provide
some historical context. Section 2.2.1 describes the the configuration of the feedback
loop, and describes the difference between this experiment, which incorporates noise
from individual photons, and earlier experiments. Section 2.3 describes the digital
signal processing and data acquisition in this experiment. Section 2.4 describes how
the photon statistics change from Poisson statistics to low-dimensional chaos as
the photon rate increases, and how different scales of observation resolve different
features of the dynamics. The next chapter discusses entropy production in this
8
system in the wider context of random number generation.
2.2 Motivation and background
In many situations, continuous variables are used to model processes which are
comprised of many individual, random occurrences. For example it is common to
model the birth and death of organisms in a population by differential equations [39].
Another example is the arrival of individual photons at a detector [40], which is often
both modeled and measured as a continuously variable intensity. It is often the case
that both large scale deterministic fluctuations and small scale noise are present.
We have constructed a model system in which small-scale shot noise from individual
photon arrivals is coupled to large scale chaotic dynamics. We describe how the
unpredictability of this system depends on the scale at which it is observed. At
small scales, we resolve individual photon arrivals, and at large scales we resolve low-
dimensional deterministic dynamics. This chapter quantifies the scale dependence
in terms of the autocorrelation function, the distribution of photon arrivals in a time
window, and the dependence of the variance of this quantity on the time window.
The next chapter describes how entropy production in this system reflects both
processes.
The experimental system we discuss in this chapter is modeled after earlier
experiments with optoelectronic feedback loops. Optoelectronic feedback loops have
become a widely-used test bed to study many concepts in nonlinear dynamics [17].
They have several features that make them appealing as a model system: the
9
timescales, dimensionality and Lyapunov spectra of their dynamics can be tuned
over a wide range, they are well-modeled by a fairly simple delay-differential equa-
tion, and they can be constructed from widely-available telecommunications compo-
nents. The wide bandwidth and separation of timescales possible in these systems
allows for breathers [19, 20], and a multistability between a large number of peri-
odic modes [41]. The high dimensionality of these systems allows them to exhibit
chimera states [38], and allows for their use in reservoir computing [42]. Opto-
electronic feedback loops have been used to demonstrate high speed chaos-based
communication over existing fiber optic telecommunications infrastructure [43]. Fi-
nally, optoelectronic oscillators have been used extensively as a model system to
study the synchronization of networks of coupled chaotic systems [18,44–46].
The complex dynamics displayed by optoelectronic feedback loops are due to
four features of their dynamics: feedback with gain, time delay, filtering, and non-
linearity. These are common features of a variety of oscillators. Feedback with gain
is a way of continually injecting energy into a system to ensure that the dynam-
ics do not decay into a trivial state. Feedback loops will typically have a finite
bandwidth, and these bandwidth limitations can often be modeled as a filter. We
introduced a filter using digital signal processing, but this is not necessary in all
situations [19]. Time delay, which we denote here as Td, is often due to finite propa-
gation time for signals. Some implementations of optoelectronic feedback loops use
long optical fibers or coaxial cables [17,19]. Here, we implement delay using digital
processing. Time delay allows for high-dimensional dynamics. The phase space of a
time-delayed system is infinite-dimensional, because, in order to determine the time
10
evolution of a system with time-delayed feedback, one must know the history of its
previous behavior over an interval of time. The initial condition is a function of
time, rather than a point in phase space. Finally, there must be some nonlinearity
present in order for the existence of chaos. Nonlinearities are usually present in
feedback loops. If the feedback gain is too low, the loop cannot oscillate, and small
perturbations will decay to a trivial fixed point. If the feedback gain is large enough,
oscillations will grow until their growth is limited by some mechanism. Such a mech-
anism necessarily “turns on” when the amplitude of oscillations is large. This kind
of amplitude dependence is not found in linear differential equations. In our opto-
electronic feedback loop, a Mach-Zehnder modulator (MZM) imposes a sinusoidal
nonlinearity.
The construction of optoelectronic feedback loops is documented in detail in
the PhD theses of Adam Cohen [47], Bhargava Ravoori [48], Caitlin Williams [49],
as well as in several papers [17,18]. In the interest of brevity, I will not go into great
detail on the aspects of this work that are well-documented by earlier publications.
Instead, I will focus on what is new in this work, namely feedback based on discrete
photon arrivals.
2.2.1 Discrete and analog dynamics
The optoelectronic feedback loop we describe in this chapter differs from the
previous experimental systems in that it incorporates a photon-counting detector
whereas other experiments incorporate analog photodetectors. Thus, the discrete,
11
photon counting
detector
variable 
attenuator
850 nm laser
electro-optic
modulator
FPGA
amplifier
DAC
time delay, Td filter
T2
T1
electrical
signal
optical 
signal
detectorlaser
electro-optic
modulator
amplifier
time delay filter
a)
b)
time
fil
te
r o
ut
pu
t
de
te
ct
or
 o
ut
pu
t
de
te
ct
or
 o
ut
pu
t
time
fil
te
r o
ut
pu
t
Figure 2.1: Schematic configuration of optoelectronic feedback loops
with an analog photodetector a), and photon-counting detector b). Qual-
itative sketches of the temporal behavior of the detector output and the
filter output are also shown. With an analog detector, the signal from
the detector is a smooth function of time. The signal from the photon-
counting detector is a series of sharp voltage pulses. The filter outputs
are roughly similar on a large scale between the two systems, but the
discreteness of the photon counts is visible.
12
stochastic individual photon arrivals are coupled with the deterministic dynamics
of the feedback loop.
Figure 2.1 shows a schematic the two types of optoelectronic feedback loops,
both with an analog detector (2.1a) and a photon-counting detector (2.1b). Light
from a fiber-coupled continuous-wave laser passes through the MZM. The MZM
modulates the optical power of the transmitted light (P ), according to the applied
voltage v, so that there is a sinusoidal relationship between the transmitted optical
power and the applied voltage: P = P0 sin2 (αv + φ). The optical power transmitted
through the modulator is then detected by a photodetector.
In the analog feedback loop shown in Figure 2.1a, the voltage signal from the
photodetector is proportional to the optical power vmod(t) ∝ P (t). This signal is
delayed and processed by a filter, and the filter output v(t) drives the modulator.
The equations of motion for a deterministic optoelectronic feedback loop can be cast
in this form [18]:
dx
dt
= Ex + βF I(t)
I(t) = sin2
[
GTx(t− Td) + φ
]
(2.1)
Here, x is the state variable of a linear, time invariant filter, matrix E and
the vectors F and G describe the characteristics of the filter, I(t) is the normalized
intensity of light transmitted through the MZM, and Td is the time delay. The
round-trip gain is denoted by β. When a photon-counting detector is used in place
of an analog photodiode, as shown in Figure 2.1b, the filter variables can be modeled
by a linear differential equation driven by discrete photon arrivals.
13
dx
dt
= Ex + βF
1
λ0
∞∑
i=1
δ(t− ti) (2.2)
When modeling this system, the photon arrivals times, {ti}, are generated by
a non-stationary Poisson point process whose rate, λ(t), depends on the state of the
filter variables.
λ(t) = λ0I(t) = λ0 sin
2
[
GTx(t− Td) + φ
]
. (2.3)
Figure 2.1 also qualitatively compares the dynamics of the filter variables as
given in Equations (2.1) and (2.2). In the deterministic system shown in Figure 2.1a,
both the time variation of the intensity and the filter output are smooth functions.
In the photon-driven system shown in Figure 2.1b, the detector produces sharp
voltage pulses whenever a photon hits it. The filter output is driven by these pulses,
the discrete nature of the photon arrivals is reflected in the sharp edges of the filter
output. Still, the shape of the filter output in Figure 2.1b is qualitatively similar
to Figure 2.1a, at least on a coarse scale. One might expect that as the incident
photon rate λ0 increases, the dynamics of these two systems become increasingly
similar.
2.2.2 Components of the feedback loop
In our experiment, the light source is a narrow-linewidth distributed-feedback
laser with a wavelength of 850 nm. The laser is operated well above threshold, and
produces a continuous wave (CW) output. The laser is coupled into a single-mode
(SM) optical fiber.
14
After passing through the modulator and attenuator, the light is coupled into
the photon counting detector. The laser wavelength of 850 nm was chosen be-
cause silicon detectors have high sensitivity in this range. We used a Geiger-mode
avalanche photodiode. The dead time of this device was about 100 ns, and the dark
count rate was about 100 counts/s. The dead time and dark counts limit the dy-
namic range of the detector. We can resolve photon rates from about 102 counts/s
to 107 counts/s.
The light from the laser is coupled into a Mach-Zendher modulator. The power
of the light exiting the modulator is related to the applied voltage by
P (t) = P0 sin
2
(
piv
2Vpi
+ φ
)
. (2.4)
For our modulator, Vpi is 1.6 V. The rate of photon arrivals at the detector is given
by:
λ(t) = λ0 sin
2
(
piv
2Vpi
+ φ
)
(2.5)
In our implementation, the filter and time delay are implemented digitally
using a field-programmable gate array (FPGA). The FPGA is an Alterra Cyclone
II, and it is incorporated into a KNJN Saxo-Q acquisition board. In addition to the
FPGA, the Saxo-Q board has two 10-bit digital-to-analog converters (DACs) with a
maximum speed of 165 MSa/s, as well as four digital-to-analog converters (ADCs)
with a maximum speed of 200 MSa/s. This board also has a USB2 interface, which
is used to configure the device, and communicate between a PC and the FPGA. We
use the USB to set the values of filter variables and transmit photon arrival time
15
data to the computer.
The amplifier is a simple non-inverting amplifier with a feedback resistance of
27 kΩ, a resistance to ground of 1 kΩ, and a fixed gain of 27. Signals with sharp
edges cause this amplifier to introduce ringing artifacts. To avoid this ringing, we
used the two-pole filter described in the next section, whose impulse response does
not have a sharp rising edge.
2.3 Digital Implementation on FPGA
The filtering in our system, as well as the time delay, are implemented digitally
on the FPGA. In order to perform these operations, the FPGA must record each
photon arrival. The photon arrival times offer the highest possible resolution of the
time evolution of the system. For this reason, the FPGA is configured to transmit
the photon arrival time data to a PC over a USB connection so that it can be stored.
It is also convenient to have the FPGA be able receive the values of parameters over
USB, so that these values can be set by computer without hard-coding them into
the verilog code. This section describes how the FPGA detects photons, updates the
filter variables, performs a time delay, and transmits photon data to the computer
on a high level. Major design choices are justified, but low-level details of the
implementation are not.
16
2.3.1 Photon Detection and Time delay
When a photon is detected by the photon-counting detector, the detector
generates a TTL pulse. We connected the output of the photon-counting detector
to an input pin of the FPGA, and the FPGA can be configured to update the value
of a register whenever the detector generates a pulse. In practice, it was helpful to
implement a deadtime in the software running on the FPGA, in order to prevent
the same photon from being registered more than once. The detected photon data
is then digitally time delayed, in order to introduce a time delay into our feedback
loop. This delay was accomplished using a shift register.
In digital signal processing, all quantities are updated in discrete time, usually
dictated by the rising or falling edge of a clock signal. A shift register is a convenient
way to implement a delay in discrete time. A shift register consists of a series of
locations in memory: Y0 . . . YNd−1, where Nd is the time delay in clock cycles. The
value at the ith location at the end of the nth clock cycle is denoted Yi[n]. On each
clock cycle, each location except for YNd−1 is updated by shifting the value to an
adjacent location in memory.
Yi[n+ 1] = Yi+1[n] for i = 0 . . . Nd − 2 (2.6)
The first location, YNd−1, is updated based on whether a photon is detected or not.
YNd−1[n+ 1] =



1 if one or more photons arrived in iteration n
0 else
(2.7)
The result is that the if system has been running for at least Nd clock cycles (one
17
delay time), we have
Y0[n] =



1 if one or more photons arrived in iteration n−Nd
0 otherwise.
(2.8)
The time-delayed photon data is represented by Y0[n], and this signal is used
to update the filter variables. I used an array of size Nd = 32768. This leads to a
maximum time delay of Td = ∆t = 1.734 ms, where ∆t is the clock cycle which is
used to update the shift register, ∆t=8/(151.1515 MHz) = 52.9270 ns.
In this experiment, two different clock signals are used for photon detection.
The clock which measures the photon arrival times has a frequency of 151.1515
MHz. The clock that is used for signal processing is 8 times slower, with a speed of
18.8939 MHz. I chose this clock speed because the dead time of the detector limits
the number of photons that can be detected in one delay time. With ∆t = 52.9270
ns, it is possible to have two photons arrive in adjacent time bins, and so the dead
time of the detector does not impose a limitation on the dynamic range of the loop.
2.3.2 Filter equations
The filter state variables are represented by 32-bit integers, and fixed point
arithmetic is used to calculate their time evolution. For readability, in this section,
I will adopt the convention that the filter variables take values between 0 and 1,
although in their digital representation, they take values between 0 and 232 − 1.
The filter state-space variables are updated according to the following equa-
18
tions.
X1[n+ 1] =
(
1−
∆t
T1
)
X1[n] + gY0[n]
X2[n+ 1] =
(
1−
∆t
T2
)
X2[n] + gY0[n]
(2.9)
where g is a factor that sets the amplitude of the filter variables, and T1 and T2
are the time scales of the filter. The DAC outputs a positive voltage signal, so that
the voltage applied to the modulator is v[n] = V0(X1[n]−X2[n]). V0 is a constant,
measured in volts, that describes the characteristics of the DAC and the amplifier.
The scale of the variables X1 and X2 is arbitrary. To put these variables in a
physical context, we note that the MZM is driven by the voltage v[n] = V0(X1[n]−
X2[n]). So, X1 and X2 should be compared to Vpi. Introducing the dimensionless
variables xi = gV0Xi pi2Vpi and β = gV0λ0
pi
2Vpi
, it is straightforward to verify that
equation (2.9) is the Euler discretization of
d
dt




x1
x2



 =




− 1T1 0
0 − 1T2








x1
x2



+ β




1
1




1
λ0
∞∑
i=1
δ(t− ti), (2.10)
which is of the same form as equation (2.2). The rate of photon arrivals is given by
λ(t) = λ0I(t) = λ0 sin
2 [x1(t− Td)− x2(t− Td) + φ] . (2.11)
2.4 Dynamics
If we take the expectation value of the random term in equation (2.10), we
obtain the deterministic model, given in equation (2.1). This model has rich dy-
namics. It is well-known that equations of the form (2.1) will go through a series
19
bifurcations into to chaos as β increases [18]. Recall that in our photon-counting
feedback loop,
β = gV0λ0
pi
2Vpi
(2.12)
So, if the quantity gλ0 is kept constant as λ0 is varied, one will see dynamics whose
noisiness varies with the values of λ0, but the analogous deterministic equation as
given by equation 2.1 will be unchanged. It is also well-known that the dimension-
ality of chaotic systems with time delay can be tuned over a wide range, and that
the dimension and number of positive Lyapunov exponents tends to increase with
the time delay [50].
In this section, we will describe the dependence of the dynamics on the incident
photon rate, λ0 with βTd held constant at 8.87. In our implementation, the time
delay is Td=1.734 ms. This is set by the memory of the FPGA and the filter clock
frequency ∆t. For the data presented here, the modulator bias is φ = pi/4, and
the filter constants T1=1.2 ms, and T2=60 µs. We vary the photon rate over a
factor of 256, from λ0Td = 12.5 (7.20 × 103 count/s) to λ0Td = 3, 200 (1.845 × 106
counts/s). As the photon rate is increased, we see a transition from Poisson shot
noise to deterministic chaos.
2.4.1 Time series, probability distributions, and autocorrelation func-
tions
Figure 2.2 shows several time series recorded with this system with increasing
photon rate. We plot Nw(t), the number of photon arrivals in the interval [t−w,t]. In
20
C
(t
’)
0
1
N
w
(t
)
0
4
8a)
C
(t
’)
0
1
N
w
(t
)
0
60
40
20
b)
C
(t
’) 1
0
N
w
(t
)
0
400
800
c)
time shift, t’ (ms)
-10 0 100 10 20 30
time t (ms) P(N) (a.u.)
C
(t
’)
0
1
I w
(t
)
0
1
0.5
d)
λ
0
T
d 
= 12.5
λ
0
T
d 
= 200
λ
0
T
d 
= 3200
deterministic simulation
Poisson
Figure 2.2: Time series, probability distributions, and autocorrelation
functions. a) - c) show experimental data, and d) shows the result of a
deterministic simulation of Equation 2.1. a) At a low photon rate of λ0Td
= 12.5, the dynamics appear Poissonian. The time series has no visible
structure, the autocorrelation function is sharply peaked at 0, and the
distribution of photon counts in a window of w = Td/4 is nearly Poisson.
b) λ0Td = 200, a slow modulation of the photon rate is evident. c) at λ0Td
= 3,200, the photon rate varies smoothly, the photon count distribution
is bimodal and much wider than a Poisson distribution with the same
mean, and the autocorrelation function shows slow oscillations. The
deterministic simulation d) shows the same features as the high photon
rate data shown in c).
21
Figure 2.2, all of the plots were generated with w = Td/4. When the incident photon
rate is λ0Td = 12.5, the photons appear to arrive at random, uncorrelated times as
in a stationary Poisson process. Increasing the incident photon rate to λ0Td = 200,
a smooth modulation of the photon rate starts to become apparent. At λ0Td =
3, 200, Nw(t) has a smooth character, and qualitatively resembles a low-dimensional
chaotic signal. We also plot the results of a deterministic simulation using equation
(2.1). This time series was smoothed with a moving average over a time window
of width w to be directly comparable with Nw(t). We plot the autocorrelation
function, C(t′) =
〈
(Nw(t)− N¯w)(Nw(t− t′)− N¯w)
〉
, normalized so that the value of
the autocorrelation function is unity at t′ = 0. As the photon rate increases from
λ0Td = 12.5, the autocorrelation function changes from a δ-like peak, characteristic
of a Poisson process, to an oscillatory function which shows correlations at long
timescales (tens of milliseconds). The autocorrelation function of the deterministic
simulation time series is in close agreement with the autocorrelation function of the
photon arrivals with λ0Td = 3, 200. Histograms of Nw(t) also show a transition from
a nearly Poisson distribution to a bimodal distribution. The peaks of the bimodal
distribution are due to turning points imposed by the maximum and minimum
values of the transmitted light through the modulator. This bimodal distribution is
seen in the deterministic simulation.
22
2.4.2 Attractor visualizations
To visualize the development of chaos with increasing photon rate, we show
Poincare´ surfaces of section in Figure 2.3. We perform a time-delay embedding of
the experimental time series Nw(t), using a time delay of ∆ = Td/4, by constructing
a list of points in two dimensional space of the form [Nw(t), Nw(t−∆)]. Because
the attractor has a dimension higher than 2, we reduce the dimensionality of the
attractor by plotting the points only when the state variables pass through a codi-
mension 1 Poincare´ surface defined by x1− x2 = pi. The embeddings show a similar
trend to the plots in Figure 2.2. We see a development of complex chaotic dynamics
from discrete photon noise as the photon rate increases. The deterministic simula-
tion is plotted for comparison, and, as in Figure 2.2, a moving average of width w
is employed so that the smoothed intensity time series, Iw(t), is directly analogous
to Nw(t). The deterministic signal in Figure 2.3d can be regarded as the infinite
photon rate limit of the photon counting system.
2.4.3 Variance
Figure 2.4 shows the dependence of the variance of Nw on the window w and
offers another indication of the transition from shot noise to deterministic chaos.
The time integral of an uncorrelated random signal executes a random walk in
which the variance grows linearly with the integration time. For this reason, we plot
Var(Nw)/w in Figure 2.4. We see distinct asymptotic growth rates of the variance
with small and large w. When w is small, the variance reflects the Poissonian nature
23
0
20
2
4
4 6
6
8
8
0
60
40
20
0
80
60 804020
800
600
400
200
0
8006004002000 10.750.50.250
1
0.75
0.5
0.25
0
Nw(t)
N
w
(t 
- ∆
)
N
w
(t 
- ∆
 )
Nw(t)Nw(t)
N
w
(t 
- ∆
 )
Iw(t)
I w
(t 
- ∆
)
Pr
ob
ab
ilit
y 
(a
.u
.)
λ0Td = 12.5 λ0Td = 200
λ0Td = 3200 deterministic simulation
a)
c)
b)
d)
Figure 2.3: Poincare´ sections. We visualize the emergence of a chaotic
attractor from Poisson noise with increasing photon rate by embed-
ding photon count time series in two dimensions with a time delay of
∆ = Td/4, and reducing the dimensionality of the dynamics by plotting
points only when the state of the system in phase space passes though a
codimension-1 surface defined by x1−x2 = pi. a) - c) show experimental
data, and d) shows the result of a deterministic simulation. These his-
tograms are constructed with a bin width of 1 photon in in a) and b), 4
photons in c), and 0.005 in d)
24
of the photon arrivals, and the growth rate of the variance has roughly constant value
of Var(Nw)/w = λ0I¯. In the limit where the counting window is much longer than
the time scale of the variations in intensity, Nw(t) can be regarded as the sum of the
photon counts in many independent identically distributed intervals, and the central
limit theorem implies that the variance will grow in proportion to w. As we increase
the photon rate from λ0Td = 12.5 to λ0Td = 3, 200, we see an increasing offset
between the two asymptotic rates of growth of the variance. The variance can be
related to the photon rate, counting window, and the unnormalized autocorrelation
function, cI(t′) =
〈
(I(t)− I¯)(I(t− t′)− I¯)
〉
[40].
Var(Nw) = w





λ0I¯ + λ
2
0
∫ w
−w
dt′
(
1−
|t′|
w
)
cI(t
′)
︸ ︷︷ ︸
Θ(w)





(2.13)
The second term in Equation (2.13) accounts for the difference between the observed
variance and the variance of a Poisson process with the same rate. The quantity
Θ(w) has units of time, and measures the correlations in I(t) introduced by the
feedback. This quantity increases from 0 to an asymptotic value Θ∞ as w increases,
accounting for the shape of the curves shown in Figure 4. cI(t′) and Θ(w) can be
calculated from time series generated by deterministic simulations. In this case, we
find Θ∞ = 150 µs. The value of Θ∞ is related to the size of the intensity fluctuations,
and the rate at which Θ(w) approaches this asymptotic value is determined by the
timescales of the correlations of I(t).
25
V
a
r(
N
w
)/
w
 (
s
-1
)
103
104
105
106
107
108
109
10-8 10-6 10-4 10-2 10-0
counting time window, w (s)
λ
0
T
d 
= 12.5
λ
0
T
d 
= 200
λ
0
T
d 
= 50
λ
0
T
d 
= 3200
λ
0
T
d 
= 800
T
d
Var(N
w
) = 1
Figure 2.4: Experimental dependence of growth rate of variance on
counting time window, w. To indicate the timescale of the deterministic
dynamics we indicate Td. We also show a line indicating Var(Nw) = 1,
which roughly separates timescales over which less than one photon ar-
rives in the counting window from timescales over which many photons
arrive. Distinct asymptotic values of Var(Nw)/w are seen in the limits
of small and large w. The offset between these two values reflects deter-
ministic correlations in the photon arrival rate and grows with increasing
photon rate.
26
Chapter 3: Entropy generation in a photon-driven feedback loop
This chapter is based on work from the following publication
Aaron Morgan Hagerstrom, Thomas Edward Murphy, and Rajarshi Roy. Harvest-
ing entropy and quantifying the transition from noise to chaos in a photon-counting
feedback loop. Proceedings of the National Academy of Sciences 112(30):9258–9263,
2015.
3.1 Overview
The main goal of this Chapter is to describe entropy quantification in the
context of the experiment described in Chapter 2. The entropy content of datasets
shown in Chapter 2 will be analyzed, and I will comment on how the entropy
estimates indicate the transition from photon counting noise to chaos described in
Chapter 2. Finally, I will discuss the relevance of these results in the general context
of random number generation.
The Chapter is organized as follows. Section 3.2 describes why physical
random number generators are important, and discusses how most researchers in
this field approach the calculation of entropy rates. I will briefly discuss some of
the recommendations in the draft NIST SP800-90B guidelines for physical entropy
sources [52]. In Section 3.3, I give some basic definitions of information-theoretic
27
quantities and describe how entropy rates are calculated. In Section 3.4, I describe
how to calculate the entropy rate of a constant-rate Poisson process. The results
in this section, while very basic, should be familiar to anyone who wants to harvest
entropy from photon arrival times. They also are a point of comparison for our ex-
perimental data. Section 3.5 shows entropy calculations applied to the deterministic
model of the photon counting feedback loop. Section 3.6 describes the application
of the entropy calculations to experimental data from the photon-counting feedback
loop. Section 3.7 discusses the relevance of these experimental results to physical
random number generation in general.
3.2 Background
Random numbers are very important for Monte-Carlo techniques [53] and for
information security applications. In many cases, random numbers are generated us-
ing deterministic algorithms, and are called “pseudorandom numbers”. The appeal
of using pseudorandom numbers is clear – computers are fast, and pseudorandom
number generators allow for random number generation without the need for ded-
icated hardware. Especially for cryptographic applications, there is a real concern
that pseudorandom numbers are not good enough. A deterministic algorithm will
have reproducible results, and there is a possibility that an attacker could exploit
this predictability. For example, in 2012, security researchers showed that an alarm-
ing fraction – two out of every thousand – of RSA public keys collected from the
Internet offer no security because they were not generated using properly-seeded ran-
28
dom numbers [54]. In order to avoid security weaknesses of this type, the random
number generator used in Linux systems is continually re-seeded using entropy [55]
harvested from events like key-strokes and mouse movements. Intel has also started
to manufacture hardware-based random number generators [56]. There is still a
great deal of interest in devices that harvest entropy from physical sources, both
because these devices can have high speeds, and more importantly, because the un-
predictability of the random numbers can be claimed to be due to a fundamentally
unpredictable physical process.
A variety of mechanisms have been proposed for physical random number
generation using quantum mechanics [13,14]. Several commercial devices, which are
marketed as “quantum random number generators” are currently being sold [15].
Some commercial devices use weak light sources and photon-counting detectors.
One scheme employs a beam splitter and two detectors, so that whenever a photon
arrives, one detector or the other detects the photon, yielding either a 0 or a 1.
Another scheme involves harvesting entropy from the random time intervals between
photon arrivals [16].
There are also a variety of techniques for random number generation based on
chaos [11, 12, 57–66]. Many of these techniques employ optical systems with a high
bandwidth and fast detectors. A major advantage of these chaos-based techniques
is that they are fast. Whereas the bit rates of photon-counting random number
generators are measured in megabits per second, the bit rates of chaos-based random
number generators are measured in gigabits or terrabits per second [12]. Some
proponents of photon-counting based random number generation claim that chaos
29
is not “fundamentally” random, but the detection of individual photons is [15]. It is
often pointed out that, on a fundamental level, the chaotic dynamics of a laser system
are continually perturbed by microscopic quantum fluctuations [65, 66]. It’s worth
noting that not all chaos-based random number generation techniques are optical,
some rely on digital electronics [58]. It’s also worth noting that chaos isn’t the
only optical method of generating random bits. Like chaos, amplified spontaneous
emission can yield bit rates in the gigabit per second range [57,67,68].
The bit rate is a crucial figure of merit for random number generators. Usually,
the designer of a random number generator wants to show that their device can be
used for practical purposes, like Monte Carlo techniques or cryptography. For these
applications, it is important that they verify that their numbers are free of bias
and correlation. “Bias” refers to streams of bits where “1” or “0” account for more
than half of the bits. “Correlation” refers to a tendency to have longer sequences
of either “0” or “1”, or more generally, a situation where the value of one bit can
be partially predicted based on the values of previous bits. Ideally, each bit is
completely statistically independent of the other bits, and is equally likely to be
either “0” or “1”. In such a case, a bit stream is said to be “full entropy” or “one
bit per symbol”. The NIST [21] and Diehard [22] tests are commonly used to verify
that a stream of bits is indeed full entropy.
Physical processes are usually not completely free of bias and correlation. In a
system based on single-photon detection, there may be bias because the beamsplitter
is not exactly 50/50, or because there are some small manufacturing differences
between the two detectors. There can also be correlations due to non-ideal behavior
30
of the detectors [69]. The probability distribution of a measured quantity in a
chaotic system is usually not completely uniform, so chaotic systems usually have
some bias. They also tend to show deterministic correlations if sampled quickly [12,
65,66]. So, in practice, almost all physical random number generators employ some
“conditioning” to reduce bias and correlation. These approaches vary in complexity.
Two of the more simple techniques are taking the logical XOR of independent bit
streams [11] and numerical differentiation [60]. Taking the idea of conditioning to
its logical extreme, it is of course possible to seed a pseudorandom number generator
with a signal measured from a physical entropy source. Indeed, some authors employ
hash functions, which will output completely full-entropy noise, regardless of how
predictable their input is [70].
3.2.1 NIST Recommendations for Entropy Sources
The NIST [21] and Diehard [22] tests were not designed as a way to assess
physical random number generators. They were designed to certify that a bit stream
is free of bias and correlation. When conditioning is used, it can be difficult to claim
that the entropy rate certified by these tests is due to the physical system, rather
than the conditioning. Given that bias and correlation can be corrected using hash
functions, NIST has drafted the new SP800-90B guidelines [52] for physical random
number generators. These standards do not require that a physical noise source be
free of bias and correlations. Instead, they require (among other things) that:
1. Although the noise source is not required to produce unbiased and independent
outputs, it shall exhibit probabilistic behavior; i.e., the output shall not be
definable by any known algorithmic rule.
31
2. The developer shall document the operation of the noise source; this documen-
tation shall include a description of how the noise source works and rationale
about why the noise provides acceptable entropy output, and should reference
relevant, existing research and literature.
Furthermore, the new standards allow for conditioning, and several specific hash
functions are approved types of conditioning. In order to justify the use of a hash
function for conditioning, one must be able to assess the entropy rate of the noise
source:
1. When an input string with m bits of assessed entropy is provided to an ap-
proved conditioning function with an n-bit output, the resulting assessed en-
tropy is uniformly distributed across the entire n-bit output. Note that if
m ≥ n, full entropy output is not necessarily provided; see item 2.
2. When an input string with 2n bits (or more) of assessed entropy is provided
to an approved conditioning function with an n-bit output, the resulting n-bit
output is considered to have full entropy.
So, to summarize, one must assess how much entropy one is entitled to take
from a noise source, and it is not as important whether the noise source has bias and
correlations. This chapter is describes one possible route towards justifying such an
assessment.
3.2.2 Entropy Generation in Physical Systems
Entropy is usually computed from discrete quantities. Physical variables like
voltage, current, and optical intensity are often either continuous, or their discretiza-
tion is so fine (individual photons and electrons) that it makes sense to think of
them as continuous. Bits, by contrast, are discrete – either 0 or 1. Measurements of
physical systems usually involve digitizing or discretizing a quantity to some finite
32
ετ
Figure 3.1: Illustration of an time trace from a physical system sampled
to resolution in time τ , and resolution in value ε. The vertical grid
lines represent sampling times, and the horizontal grid lines represent
bin edges.
resolution, ε. Measurements of a physical system are usually discrete in time as
well. For example, an oscilloscope records samples at discrete times, with a time
interval between measurements of τ . These three parameters ε and τ are often un-
avoidable. In some cases, on might need to consider a third variable, w. In the case
of the photon-counting feedback loop described in Chapter 2, w is the time interval
over which we count photons. Rate fluctuations much faster than 1/w will not be
resolved.
33
In the context of statistical mechanics, information entropy is a very important
concept [71–75]. There is a rich literature on the dependence of the entropy rate,
h, of signals harvested from physical systems on the variables ε and τ [10, 76–80].
In these contexts, the entropy rate is often called the (ε, τ) entropy, and is denoted
h(ε, τ). It has been shown in Brownian motion and Johnson noise in RC circuits that
statistical measurements of this entropy can be related to heat dissipation [81, 82].
There have also been experimental investigations of entropy rates in convection
experiments [77,80]. The (ε, τ) entropy will have qualitatively different dependence
on ε and τ depending on the physical origin of the unpredictability [10,76,80].
3.3 Entropy and Entropy Rate
This section discusses a few basic concepts from information theory that will be
employed in the rest of the chapter, with an emphasis on calculating the entropy rate
of a dynamical system. Most of these concepts are described in standard textbooks,
like [83]. For a discussion of entropy in the context of dynamics, see [10, 76, 80, 84],
especially [10, 80].
This section requires some notation. This is all standard notation, but for
clarity I’m making a table of definitions (3.1).
The Shannon entropy of a random variable X is defined as
H(X) = −
∑
x∈X
p(x) log2 p(x). (3.1)
34
Table 3.1: some notation
Symbol Meaning
X a generic random variable
x a specific value that X might take, like 4
X the set of possible values x might take
p(x) the probability that X takes the specific value x
Xi a random variable, the value of X at time i
xi the value of Xi in a particular realization
x(t) the value of X at a specific time in
a specific realization of the physical process
p(xi, xj) the joint probability that Xi takes
the value xi and Xj takes the value xj
p(xi|xj) the conditional probability that Xi
takes the value xi given that Xj takes the value xj
The entropy of this joint probability distribution can be calculated in the same way:
H(X1, X2, . . . , Xd) = −
∑
x1∈X
∑
x2∈X
· · ·
∑
xd∈X
p(x1, x2, . . . , xd) log2 p(x1, x2, . . . , xd).
(3.2)
Another important concept is conditional entropy. This is used to quantify to
what extent one variable, X1, can be predicted by sampling a different variable, X2.
Conditional entropy is defined as the expected value of the entropy of X1 given that
X2 has been sampled and is known to have the value x2:
H(X1|X2) = −
∑
x2∈X
p(x2)H(X1|X2 = x2)
= −
∑
x1∈X
∑
x2∈X
p(x1, x2) log2 p(x1|x2)
(3.3)
If X1 can be completely predicted by sampling X2, then H(X1|X2) = 0. On the
other hand, if X1 and X2 are independent, H(X1|X2) = H(X1). The conditional
35
entropy can be related to the entropy of a joint probability.
H(X1|X2) = −
∑
x1∈X
∑
x2∈X
p(x1, x2) log2 p(x1|x2)
= −
∑
x1∈X
∑
x2∈X
p(x1, x2) log2
p(x1, x2)
p(x2)
= H(X1, X2)−H(X2)
(3.4)
This is called the “chain rule”.
The entropy rate of a stochastic process can be measured either in bits per
sample, or bits per second. For consistency with the rest of the chapter, I will write
the formulas for the entropy rate in bits per second, by adding a factor of 1/τ , where
τ is the time interval between samples. There are two definitions of the entropy rate
h. These definitions will yield the same result for stationary processes [83], but, in
practice, one is more easy to calculate.
The first definition is
h1(X) = lim
d→∞
1
τd
H(X1, X2, . . . , Xd). (3.5)
In Equation 3.5, one considers patterns of increasing length d, and measures how the
entropy of the set of patterns of length d changes with d. If d is large enough, one
finds that the entropy of the set of patterns (pattern entropy) grows proportionally
to d. This linear growth is shown in Figure 3.2.
The second definition of entropy rate is based on the extent to which the
current sample, xi, can be predicted based on the history of previous samples.
h2(X) = lim
d→∞
1
τ
H(Xd|Xd1 , . . . , X1). (3.6)
36
1.5
1.0
0.5h 
(b
its
/s
am
pl
e)
d
0 5 10 15 20
d
0 5 10 15 20
10
8
6
4
2
0
h2
h1
H
(X
1 
, X
2 
, .
.. 
, X
d 
) (
bi
ts
)
0.0
a.
b.
hτ
Figure 3.2: Pattern entropy and entropy rate. a) An example of the
growth of the pattern entropy, Hd = H(X1, X2, . . . , Xd), as a function
of d. The asymptotic linear growth is indicated by a dashed line. Al-
though this figure is intended as an illustration of concepts of pattern en-
tropy and entropy rate, experimental data are shown here. The pattern
entropy rate was calculated by applying the Cohen-Procaccia method
(Equation 3.9) to the time series Nw(t), from the photon-counting feed-
back loop, with λ0Td = 3200, τ = w = Td/4, and ε = 400 photons.
Note that there are two linear scaling regions of Hd. Small embedding
dimensions, or pattern lengths, are insufficient to predict the dynamics,
leading to a fast increase of Hd with d. b) Approximations of the entropy
rate based on Equations 3.5 and 3.8. Both h1 and h2 approach the same
limit, but h2 approaches much more quickly with increasing embedding
dimension.
37
As mentioned previously h1 = h2 for stationary systems. Equation (3.6) be re-
written in another form that is more convenient experimentally. By the same logic
as Equation (3.4), we can write
H(Xd|Xd1 , . . . , X1) = H(Xd, Xd−1, . . . , X1)−H(Xd−1, Xd−2, . . . , X1), (3.7)
which leads to
h2(X) = lim
d→∞
1
τ
[H(Xd, Xd−1, . . . , X1)−H(Xd−1, Xd−2, . . . , X1)] . (3.8)
Like Equation 3.5, this formula for the entropy rate indicates a linear growth of the
pattern entropy with d. Whereas Equation 3.5 calculates the average entropy per
symbol in a large pattern, Equation 3.8 calculates the slope of the growth of the
pattern entropy with d. It can be argued that h2 will give a better approximation
of the entropy rate for finite d [85], and Figure 3.2 suggests that this is the case.
3.3.1 Entropy Calculations
The first step to computing the (ε, τ) entropy rate of an experimental signal is
to generate a list of points in d-dimensional space using time delay embedding with
a delay of τ . These vectors can be regarded as samples of a d-dimensional prob-
ability distribution over phase space. The entropy of this probability distribution,
Hd = H(X1, X2, . . . , Xd). In principle Hd can be calculated by building a histogram
with boxes of width ε and applying Shannon’s formula, Equation 3.1 [83]. In prac-
tice, direct application of this approach requires a very large amount of data when
the embedding dimension is large. Cohen and Procaccia proposed a more efficient
38
algorithm to estimate the pattern entropy [85] in the context of estimating metric
entropy from experimental data. First, one randomly selects a small number M of
reference points from the time series. M is usually small compared to the length of
the time series. In our case we observed only small differences between M = 5000
and M = 20000. For each reference point i, one computes ni(ε), the fraction of
the points in the dataset within a box of width ε centered on the reference point.
The only difference between a direct calculation of the Shannon entropy, and the
Cohen-Procaccia procedure is that in a direct calculation, a rectangular array of
bins is used, rather than a set of bins centered on random points chosen from the
data set. In searching for neighbors for the ith reference point, we exclude points
within a time window of τ of that point, as suggested by Theiler [86]. The pattern
entropy is then estimated by
Hd(ε) = −
1
M
M∑
i=1
log2 ni(ε). (3.9)
It is a general feature of unpredictable signals that Hd grows linearly with d in the
limit that d is large, and the entropy rate is the slope of this linear increase, as indi-
cated in Equation 3.8. Cohen and Procaccia argued that Equation 3.8 approaches
the correct value of the entropy rate at a lower dimension than Equation 3.5. Figure
3.2 shows an anecdotal example where this is the case.
We will see that the calculation of Hd from Equation 3.9 requires a large
amount of data. Some authors have suggested using the correlation entropy, H(2)d ,
in place of Hd. [10, 76,77].
H(2)d (ε) = − log2
1
M
(
M∑
i=1
ni(ε)
)
(3.10)
39
3.4 Entropy rate of a constant rate Poisson process
In this section, I give a few basic results about the entropy rate of a Poisson
process. This is relevant to the discussion of a photon-counting optoelectronic feed-
back loop because the arrival of photons generated by a laser at a photon-counting
detector is described by a Poisson process. The results given in this section are far
from new, but will provide context for later results, and it can be helpful to show
how the methods described here can be applied in a simple case.
A Poisson process describes events (like photon arrivals) that occur at random
times. The probability per unit time for an event to occur is constant, and the
probability that an event occurs at time t is independent of whether or not an event
occurs at time t′.
Because the events in different points in time are statistically independent, it’s
very easy to evaluate the entropy rate. I will assume that we have divided time into
small bins of width τ , and that it is possible to sample system much faster than the
event rate, so τλ0  1. I will consider the case where τλ0  1 later in this section.
The probability that one or more events occurs in a time bin is 1 − eλ0τ ≈ λ0τ . I
will now consider a random variable Yi, which takes a value 1 if one or more photons
arrive in bin i, and 0 if no photons arrive. The probability mass function of Yi is
p(Yi = 0) = 1− p0
p(Yi = 1) = p0
(3.11)
The entropy rate can be trivially evaluated from either Equation 3.5 or Equa-
40
tion 3.6.
h =
−1
τ
[(1− p0) log2 (1− p0) + p0 log2 (p0)] (3.12)
In photon counting experiments, it is possible to time photon arrivals to a high
degree of precision. We used a clock with a speed of 1/τ=151.1515 MHz, and state-
of-the-art detectors can tag photon arrivals to tens of ps. On the other hand, the
maximum count rates resolvable by detectors are on the order of 107 counts/second.
So, we are interested in the limit that τ → 0. Recalling that p0 ≈ λ0τ , and expanding
to lowest order, we have
h = λ0
(
λ0
ln 2
− log2 λ0τ
)
(3.13)
There are a few comments to make here. The entropy rate grows without
bound as τ → 0, but the growth is logarithmic. This presents a limitation to the
scalability of random number generation using photon-counting detectors. To put
in real numbers, the experimental system described in Chapter 2 has a temporal
resolution of τ=6.62 ns. If λ0 = 2 × 106 counts/s, the entropy rate is h = 15.4
Mbits/s. If the timing resolution could be increased to τ=66.2 ps, which is on the
order of the state of the art, then we would have h = 28.6 Mbits/s. Improving the
timing resolution of a detector by a factor of 100 will not increase the bit rate by a
factor of 100 but instead will yield a relatively small increase of about 6.6 bits per
photon. In contrast, when either a chaotic system or amplified spontaneous emission
is used to generate random numbers, the entropy rate can scale more favorably with
τ . As discussed in Section 3.5, the entropy rate of a chaotic system is determined by
41
its positive Lyapunov exponents, which can be on the order of hundreds of gigabits
per second [12,57]. These fast dynamical timescales allow one to take advantage of
fast digitizers to harvest entropy at high rates.
Another point to make here is that bias and entropy rate are distinct concepts.
As τ → 0, the distribution of Yi is increasingly biased in favor of 0, and yet the
entropy rate grows.
3.4.1 Entropy rate from inter-arrival intervals
Another way to harvest entropy from a Poisson process would be to use the
inter-arrival intervals. Time will still be discretized in bins of width τ , but now
the random variable we will consider is Ti, the time interval between the (i − 1)th
photon arrival and the ith photon arrival.Since time is divided into bins of width τ ,
Ti takes on discrete values: τ , 2τ , and so on.
The variables Ti are independent and identically distributed, and this distri-
bution is given by a geometric distribution. This is a direct consequence of the
statistical independence of different time bins, and the assumption that the rate is
constant.
p(ti = kτ) = p0(1− p0)
k for k=1,2,...,∞ (3.14)
From 3.5, we have
h =
−1
t¯
∞∑
k=1
p0(1− p0) log2 p0(1− p0) (3.15)
Where t¯ = τ/p0 is the average of ti. This sum can be evaluated in closed form, and
42
we have
h =
−1
τ
[(1− p0) log2 (1− p0) + p0 log2 (p0)] , (3.16)
which is identical to Equation 3.12. So, measuring inter-arrival times will yield
exactly the same entropy rate as measuring Yi.
3.4.2 Can intensity modulation increase entropy rates?
So far, I have only discussed constant-rate Poisson processes. In the previous
chapter, I described a system in which feedback is used to modulate the intensity
of a weak light source, so that variations in the photon arrival rate were seen at
the detector. It is natural to ask if these rate fluctuations increase the entropy rate
of the signal compared to a constant-rate process. It turns out that the Poisson
process has the highest possible entropy of any process consisting of a sequence of
discrete events with a given average inter-event interval.
This can be shown using Lagrange multipliers. We will assume that we are
harvesting entropy from the intervals between events, Ti, as in the last section. The
average inter-event interval will be denoted t¯. It may be the case that when rate fluc-
tuations are present, the inter-event intervals are no longer independent or identically
distributed. I will assume that the process which determines the inter-arrival inter-
vals is stationary, but individual photon arrivals are not necessarily independent.
So p(ti = k) = p(tj = k), but p(ti|ti−1, ti−2, . . . , ti−d) 6= p(ti|tj−1, tj−2, . . . , tj−d).
It is well-known [83] that for any random variables X1, X2, . . . Xd,
H(X1, X2, . . . Xd) ≤ H(X1) +H(X2) + · · ·+H(Xd) (3.17)
43
From Equation 3.5, we then have
h ≤
−1
t¯
∞∑
k=1
p(k) log2 p(k) (3.18)
I will now maximize the right-hand-side of 3.18 subject to the constraints that the
distribution has a well-defined mean, is defined for integers k greater than 1, and is
normalized. Introducing the Lagrange multipliers α and β, the variational equations
are:
∂
∂p(k)
[
α
(
t¯− τ
∞∑
k=1
kp(k)
)
+ β
(
1−
∞∑
k=1
p(k)
)
+
∞∑
k=1
p(k) log2 p(k)
]
= 0 (3.19)
This can be solved for pk.
p(k) = 2−ατk−1−β (3.20)
This has the form of a geometric distribution. The values of α and β are set by
the constraint equations. The result will be Equation 3.14 with p0 = 1/t¯. The
conclusion is that, if entropy is harvested from the time intervals between a series of
events with a well-defined average inter-event interval, a Poisson process will yield
the highest possible entropy rate.
3.4.3 Counting windows and ε
In Chapter 2, I used the variable Nw(t), the number of photon counts in the
interval [t−w, t], to analyze the behavior of the photon-counting feedback loop. One
advantage of this variable is that it shows the chaotic character of the dynamics,
at least as far as having similar time traces and autocorrelation functions. Another
example is that it can be directly compared to a Poisson process with a constant
44
rate. In a Poisson process, the distribution number of events, N , to arrive in an
interval [t− w, t] is
p(N) =
N¯N
N !
e−N¯ . (3.21)
In a constant-rate process, the expected number of events is N¯ = wλ0
In a constant-rate process, successive samples are uncorrelated if τ ≥ w. The
entropy rate is then straightforward to calculate from Equations 3.1 and 3.5. There
is not a closed-form solution, but in a certain limit, it is possible to make a crude
approximation that works fairly well. In the limit that λ0w  0, the Poisson
distribution (Equation 3.21) approaches the Gaussian distribution with a mean and
a variance of λ0w. If
√
λ0w  ε, the sum in Equation 3.1 can be approximated by
an integral. In this case, we can approximate the entropy rate of a Poisson process
with a counting window as
h =
1
τ
log2
(√
2piewλ0

)
. (3.22)
Note that this entropy’s dependence on  is − log2 . The logarithmic dependence
on  is common in noisy systems [76, 80, 84]. Figure 3.3 shows the comparison of
the entropy rate calculated from Equation 3.22 to the entropy rate calculated from
experimental photon count data.
3.5 Entropy rate of chaotic systems
In chaotic systems, unpredictability is due to the sensitive dependence on ini-
tial conditions. Because small perturbations grow exponentially in time, chaotic
45
200
3200
λ
0
T
d 
=12.5
50
800
Gaussian 
approximation
0
2000
4000
6000
8000
h
(ε
,τ
) 
(b
it
s
/s
)
ε
4 16 64 2561
Figure 3.3: Entropy rate calculations for constant-rate Poisson processes
with λ0Td = 3200, 800, 200, 50, and 12.5, with Td = 1.734 ms. These
data were recorded using the experimental system described in Chapter
2, but with the feedback disconnected. The solid lines were calculated
using the Cohen-Procaccia algorithm, applied to Nw(t) with w = Td/4.
Equation 3.9. The dashed lines were calculated using Equation 3.22.
46
εε
a.
b.
d=5d=10
15
20
25
h(
ε,
τ)
 (b
its
/s
)
h(
ε,
τ)
 (b
its
/s
)
750
500
250
0
1/8 1/2 2
1/8 1/2 2
750
500
250
0
hks
hks
25
20
15
d=10 d=5
Figure 3.4: Dependence of h(ε, τ) on ε for a) τ = Td/4 and b) τ = 3Td/4.
In this figure, ε has the same units as the normalized intensity I(t), which
varies between 0 and 1. As one would expect, the dependence of h on ε
is flat, as long as ε is small compared to the size of the attractor. For
both τ = Td/4 and τ = 3Td/4, the level of this plateau slowly decreases
with increasing embedding dimension d, and appears to approach hks.
47
systems generate information. The growth of uncertainty is quantified by the Lya-
punov exponents µi, and in particular the largest exponent, µ1. Positive Lyapunov
exponents and entropy rate both quantify unpredictability and there is a close rela-
tionship between these two quantities. One would expect that if a chaotic system is
sampled infrequently (τµ1  1), successive samples will be uncorrelated because of
the growth of uncertainty between measurements. On the other hand, if the inter-
val between successive samples is small (τµ1  1), one expects strong correlations
between adjacent samples and a reduced entropy per sample. Experimental and
theoretical work using semiconductor lasers has shown that these considerations are
crucial to physical random number generation using chaotic dynamics [64–66]. In
the limit that τ, ε → 0, h(ε, τ) will approach a finite value, the Kolmogorov-Sinai
(or metric) entropy, hks [10, 87–89]. The metric entropy is related to the Lyapunov
exponents, µi, by
hks =
1
log(2)
∑
µi>0
µi. (3.23)
We calculated the spectrum of Lyapunov exponents from Equation 2.1, with all of
the parameters given in Chapter 2 [50, 90]. There is only one positive Lyapunov
exponent with a value of µ1/ log(2) = 345 bits/s. The Kaplan-Yorke dimension [91]
calculated from the Lyapunov spectrum is 3.56.
The (ε, τ) entropy will have qualitatively different behavior as ε→ 0 depending
on the physical origin of unpredictability [10, 80]. In chaotic systems, the entropy
rate does not depend on either the sampling rate or the sampling resolution. This
property of chaotic systems imposes a theoretical limitation on physical random
48
number generation. Increasing the speed and resolution of a measurement device
cannot in principle increase the entropy that can be harvested from a deterministic
chaotic system beyond hks. In contrast to deterministic systems, the, entropy rate
of stochastic signals diverges like − log(ε) for finite τ [76, 80].
Figure 3.4 shows the dependence of h(ε, τ) on ε for τ = Td/4 and τ = (3/4)Td.
In this figure, ε has the same units as the normalized intensity I(t), which varies
between 0 and 1. As one would expect, h in independent of ε as long as ε is small
compared to the size of the attractor. With both τ = Td/4 and τ = 3Td/4, the
level of this plateau slowly decreases with increasing embedding dimension d, and
appears to approach hks.
Figure 3.4 shows that the convergence of the entropy estimate with increasing d
is rather slow. The comparison of Figure 3.4 a) and 3.4 b) suggests that τ = (1/4)Td
is too small to resolve the expected entropy rate, hks. Care should be taken when
applying this method. The slow convergence of the entropy estimate could lead to
too-large estimates of the entropy. Comparing the curves in Figure 3.4 a) for d = 10
and d = 15, one might imagine that the entropy rate estimate given by these curves
is close to its asymptotic value for high d, but such an estimate would be off by a
factor of two.
3.6 Entropy generation in a photon-driven feedback loop
To analyze the photon count data from the experiment described in Chapter
2, we computed Nw(t), the number of photon counts in the interval [t − w, t]. We
49
εa.
200
3200
λ0Td =12.5
50
800
Poisson
ε
b.
4 16 64 256 1024
4 16 64 256 1024
250
0
500
750
1000
1250
250
0
500
750
1000
1250
hks
hks
h(
ε,
τ)
 (b
its
/s
)
h(
ε,
τ)
 (b
its
/s
)
Figure 3.5: Entropy rate of experimental time series. The entropy is
calculated from Nw(t), and ε is measured in photons. τ = Td/4 is
shown in a), and τ = (3/4)Td is shown in b). Embedding dimensions
of d=5,10,15,20, and 25 are shown. At low photon rates we see a diver-
gence characteristic of a Poisson processes. As the photon rate increases,
the dependence of the entropy on ε becomes progressively flatter, and
approaches hks. Across photon rates, we see a divergence of the entropy
rate for small ε. The Poisson curves were calculated using Equation 3.22
with a rate of λ0.
50
also chose a rather large value of τ , on the order of Td. These two choices allowed
us to detect deterministic chaos in the experimental time series.
Figure 3.5 shows the entropy per unit time calculated from the experiment with
both τ = (1/4)Td and τ = (3/4)Td. We see that as the photon rate increases, the
dependence of the entropy on ε becomes progressively flatter at high ε. Furthermore,
in the region that this flattening is present, the value of h(ε, τ) is close to hks.
The flattening of h(ε, τ) at high photon rate is another indication that this system
behaves more deterministically in the high photon-rate regime. At all photon rates,
we see h(ε, τ) sharply increases as ε decreases, which is due to the shot noise inherent
in the system. It is natural to compare the entropy rates we observe to a constant-
rate Poisson process with the same average rate. In Figure 3.5 a), we see the entropy
rate approaching the rate of the Poisson process given by the dashed lines for small
ε. In Figure 3.5 b), it is harder to infer the asymptotic behavior for small ε, because
we were unable to resolve the value of h with sufficiently small ε.
It’s also worth noting that in comparing Figure 3.4, and 3.5, the difference
between τ = (1/4)Td and τ = (3/4)Td is more pronounced in simulation than
experiment. The entropy rate also converges faster as d increases in experiment
than in simulation.
3.7 Comments
At the outset of this project, we had a rather simple question: does discrete
noise, for example arising from individual photon arrivals, increase the unpredictabil-
51
ity of a chaotic system? Another, related question, is whether chaos “amplifies” the
discrete noise. As I describe in Chapter 2, we were able to give some insight into
these questions by examining the qualitative appearance of the time series and
Poincare´ sections (Figure 2.3), the variance of Nw(t) (Figure 2.4) and the autocor-
relation functions (Figure 2.2). The autocorrelation function was sharply peaked
at low photon rates, indicating a lack of temporal correlations. As the photon rate
increased, the autocorrelation function developed a complex, oscillatory character
with correlations lasting for tens of ms, similar to the deterministic system. The
qualitative appearance of both the time series and Poincare´ sections confirmed that,
as the photon rate λ0 increases, deterministic chaos develops from noise. So, in a
sense, The lower photon rates are noisier.
We then set out to show these trends by quantifying entropy production. It
quickly became apparent that we could not avoid considering three parameters: w
(the counting window), ε (resolution in value), and τ (time resolution) if we wanted
to understand the roles of chaos and shot noise. If we employ the smallest value
of τ possible, and count individual photon arrivals, we are left with the straight-
forward result of Subsection 3.4.2, that the constant rate Poisson process has the
highest achievable entropy. Furthermore, the entropy scales with the photon rate in
a straightforward manner, given in Equation 3.13. We are then left with an unsat-
isfying conclusion: a constant-intensity light source with a high photon rate gives
the largest entropy rate. This analysis is unsatisfying because it fails to capture the
unpredictability in the chaotic wave forms shown in Figure 2.2. It also fails to cap-
ture the comparatively noisy character of the time traces for lower photon rates. In
52
order to resolve these features of the dynamics, we had to introduce coarse graining
in both time (τ and w) and value, ε. The choice of (τ and w) allow for us to resolve
the slow intensity fluctuations and the associated uncertainty. The dependence of
the entropy on  showed that at small scales, Poisson shot noise dominates, and at
large scales, the chaotic dynamics can be resolved.
3.7.1 Broader context
Two prevalent methods have attracted attention for optical random number
generation: those based on single photon detection from strongly attenuated light
sources [13, 14], and those based on digitized high-speed fluctuations from chaotic
lasers [11]. In the former case, the entropy is claimed to originate entirely from
quantum mechanical uncertainty, yet in practice these methods are also subject to
unpredictable drift and environmental variations. In the latter case, the entropy is
attributed to the dynamical unpredictability of chaos, but the unavoidable presence
of spontaneous emission is thought to play a role in seeding these macroscopic fluctu-
ations [64,66]. The system presented here is unprecedented in that it can approach
macroscopic chaos from the single photon limit, thereby revealing the transition
from noise to chaos. Moreover, the analysis offers a unified measure of entropy that
captures both behaviors, and clarifies the relationship between sampling frequency,
measurement resolution and entropy rate.
The designer of a physical random number generator must choose the sam-
pling rate and resolution that they will use to collect numbers from a physical
53
system. These decisions will impact the entropy rate. Heuristically, finer discretiza-
tion (smaller ε) and more frequent sampling (smaller τ) lead to higher entropy rates,
but without the methods presented here it is difficult to assess the dependence of
the entropy rate on these parameters in any given system. The statistical tests that
are usually used to evaluate physical random number generation [21, 22] were not
designed to answer these questions, but rather to certify that a stream of bits is free
of bias and correlation. If a random number generator employs post processing (as
most do), existing statistical tests applied to the output binary sequence provide no
insight into whether the entropy originates from the physical process or the post-
processing algorithm employed. We offer a new approach to determine the entropy
of a physical process that clarifies the origin and nature of uncertainty and informs
the choice of sampling rate and measurement resolution.
54
Chapter 4: Cluster Synchronization and Symmetries in Networks
This chapter is based on work from the following publications
Louis M. Pecora, Francesco Sorrentino, Aaron M. Hagerstrom, Thomas E. Murphy,
and Rajarshi Roy. Cluster synchronization and isolated desynchronization in com-
plex networks with symmetries. Nature Communications, 5, 2014.
Francesco Sorrentino, Louis M. Pecora, Aaron M. Hagerstrom, Thomas E. Murphy,
and Rajarshi Roy. Cluster Synchronization in Laplacian Coupled Systems: Beyond
symmetries. Submitted to Physical Review Letters (arXiv preprint arXiv:1507.04381),
2015
4.1 Overview
This chapter focuses on synchronization patterns that can emerge in networks
of coupled dynamical systems. The symmetries of a network can be used both to
predict the synchronization patterns, and to analyze their stability.
This chapter is organized as follows. Section 4.2 provides background and
motivation, and introduces the mathematical formulation of the problem of clus-
ter synchronization. Section 4.3 describes the relationship between the symmetries
of a network and the possible cluster synchronization patterns. Section 4.4 de-
scribes, very briefly, the linear stability analysis of global synchronization. Section
4.5 describes the linear stability analysis of cluster synchronization. Section 4.6
55
describes the a coordinate transformation that is necessary for this analysis. Sec-
tion 4.7 describes the experimental system that I used to obtain the experimental
demonstration of cluster synchronization shown in Section 4.8. Section 4.9 briefly
describes the generalization of the symmetry-based approach to a much larger class
of synchronization patterns.
4.2 Background
Synchronization in complex networks is essential to the proper functioning of a
wide variety of natural and engineered systems, ranging from electric power grids to
neural networks [92]. Global synchronization, in which all nodes evolve in unison, is a
well-studied effect, the conditions for which are often related to the network structure
through the master stability function [27]. Equally important, and perhaps more
commonplace, is partial or cluster synchronization (CS), in which patterns or sets
of synchronized elements emerge [26]. Recent work on cluster synchronization has
been restricted to networks where the synchronization pattern is induced either by
tailoring the network geometry or by the intentional introduction of heterogeneity in
the time delays or node dynamics [46,93–99]. These anecdotal studies illustrate the
interesting types of cluster synchronization that can occur, and suggest a broader
relationship between the network structure and synchronization patterns. Recent
studies have begun to draw a connection between network symmetry and cluster
synchronization, although all have considered simple networks where the symmetries
are apparent by inspection [23–25]. More in-depth studies have been done involving
56
bifurcation phenomena and synchronization in ring and point symmetry networks
[100, 101]. Finally, a form of cluster synchronization can occur in situations where
input and output couplings match in a cluster, but there is no symmetry [101,102].
Here we address the more common case where the intrinsic network symmetries
are neither intentionally produced nor easily discerned. We present a comprehensive
treatment of cluster synchronization, which uses the tools of computational group
theory to reveal the hidden symmetries of networks and predict the patterns of
synchronization that can arise. We use irreducible group representations to find
a block diagonalization of the variational equations that can predict the stability
of the clusters. We further establish and observe a generic symmetry breaking
bifurcation termed isolated desynchronization, in which one or more clusters lose
synchrony while the remaining clusters stay synchronized. The analytical results
are confirmed through experimental measurements in a spatiotemporal electro-optic
network. By statistically analyzing the symmetries of several types of networks, as
well as electric power distribution networks, we argue that symmetries, clusters,
and isolated desynchronization can occur in many types of complex networks and
are often found concurrently. Throughout we use the abbreviation ID for isolated
desynchronization.
4.2.1 Coupling – adjacency matrix and Laplacian matrix
Two types of matrices are very commonly used to describe networks [103]: the
adjacency matrix and Laplacian matrix. The adjacency matrix is denoted by A.
57
The element Aij is 1 if nodes i is connected to node j, and 0 otherwise. Here, we
only discuss the case of bidirectional coupling, where Aij = Aji, but there are many
cases where one might want to consider a network in which there are unidirectional
couplings. It’s also common to have weighted connections. In this case, the value
of Aij is not necessarily 0 or 1, and specifies the strength of the connection between
i and j. Most of the formalism discussed in this chapter will probably carry over
to both weighted and unidirectional connections, but this has not been explored in
detail. Typically, the diagonal entries of A are set to 0.
The Laplacian matrix is the same as the adjacency matrix, except for the
values of the diagonal elements. In the Laplacian matrix, normally denoted L, the
diagonal entries are given by the opposite of the row sums.
Lij = Aij for i 6= j
Lii = −
∑
j 6=i
Aij
(4.1)
As a consequence of this definition,
∑
j Lij = 0. In the next section, I will discuss
how this allows for global synchronization.
Except where otherwise noted, this chapter discusses adjacency matrix cou-
pling where all of the links are bidirectional and of equal weight.
4.2.2 General model of networks of coupled identical systems
A set of general dynamical equations to describe a network of N coupled
identical oscillators are
x˙i(t) = F(xi(t)) + σ
∑
j
AijH(xj), i = 1, ..., N, (4.2)
58
where xi is the n-dimensional state vector of the ith oscillator, F describes the
dynamics of each oscillator, A is a coupling matrix that describes the connectivity
of the network, σ is the overall coupling strength, and H is the output function
of each oscillator. Eq. 4.2 or its equivalent forms provide the dynamics for many
networks of oscillators including all those in Refs. [23–27, 92, 93, 95–97, 99]. This
includes some cases of time delays in the coupling functions. As noted in Ref. [92]
it is only necessary for the form of the equations of motion or, more importantly,
the variational equations to have the form of Eq. 4.2 near the synchronization
manifolds. The form of Eq. (4.2) also applies to discrete time systems or more
general coupling schemes [104]. And the same form emerges in the more general
case for the variational equations where the vector field and the coupling combine
into one function, viz. F(xi(t), {xj(t)}), where the argument in brackets is the
input of all nodes connected to node i, so long as the nodes are treated as having
the same basic dynamics. See Ref. [104] for more explanation. Generalizations of
Eq. (4.2) have been studied in Refs. [46,94,98,105]. The types of natural and man-
made systems which can be modeled by equations of the same form as Eq. (4.2) is
large [103,106].
The same analysis as presented here will apply to more general dynamics.
Here, for simplicity, we take all nodes to have identical dynamics and be bidirec-
tionally coupled to other nodes in the network by couplings of the same weight, i.e.
Aij is taken to be the (symmetric) adjacency matrix of 1’s and 0’s with the factor
of σ controlling the weight of the couplings. The vector field F can contain self
feedback terms. Here we take the self feedback to be identical for all nodes.
59
4.3 Symmetries and cluster synchronization
Symmetries are relabelings of the nodes of a network that do not change the
equations of motion. Figure 4.1 shows a network with a few symmetries that are
apparent by visual inspection. For example, the network can be reflected about
the three different axes indicated. It can also be rotated about the central node,
either clockwise or counter-clockwise. In total, this network has six symmetries,
including the trivial symmetry (the identity). The symmetries of the network form
a (mathematical) group G. Each symmetry g of the group can be described by
a permutation matrix Rg that re-orders the nodes. These permutation matrices
commute with the adjacency matrix, i.e. A = RTg ARg for each Rg. As a result, the
dynamical equations are unchanged by these operations.
The set of symmetries (or automorphisms) [100,107] of a network can be quite
large, even for small networks, but they can be calculated from A using widely avail-
able discrete algebra routines [108, 109]. Figure 4.2a shows three graphs generated
by randomly removing 6 edges from an otherwise fully connected 11-node network.
Although the graphs appear similar and exhibit no obvious symmetries, the first
instance has no symmetries (other than the identity permutation), while the others
have 32 and 5,760 symmetries, respectively. So for even a moderate number of nodes
(11) finding the symmetries can become impossible by inspection.
Once the symmetries are identified, the nodes of the network can be parti-
tioned into M clusters by finding the “orbits” of the symmetry group: the disjoint
sets of nodes that, when all of the symmetry operations are applied, permute among
60
1
3
2
4
1
2
4
3
1
4
3
2
1
4
2
3
1
2
3
4
1
3
4
2
1 0
0 0
0 0
1 0
0 0
0 1
0 1
0 0
1 0
0 0
0 0
0 1
0 1
0 0
0 0
1 0
1 0
0 1
0 0
0 0
0 0
0 0
1 0
0 1
1 0
0 1
0 0
0 0
0 0
0 0
0 1
1 0
1 0
0 0
0 0
1 0
0 1
0 0
0 0
0 1
1 0
0 0
0 0
0 1
0 0
0 1
1 0
0 0
Figure 4.1: Symmetries and clusters in a simple network. All of the
symmetry operations of this graph are shown, as well as the permuta-
tion matrices that represent them. Note that the outer nodes, 2,3, and
4 can change position under the operations, but node 1 cannot be per-
muted with 2,3, or 4. It is possible for nodes 2,3, and 4 to synchronize
independently of 1. A solution of this form is allowed because if the
nodes synchronize in this way, the equations of motion for nodes 2,3,and
4 will be identical.
61
one another [100]. Because Eq. (4.2) is essentially unchanged by the permutations
the dynamics of the nodes in each cluster can be equal, which is exact synchroniza-
tion. Hence, there are M synchronized motions {s1, ..., sM}, one for each cluster. In
Fig. 4.2a, the nodes have been colored to show the clusters. For the first example,
which has no symmetries, the network divides into M = N trivial clusters with one
node in each. The other instances have 5 and 3 clusters, respectively. Once the clus-
ters are identified, Eq. (4.2) can be linearized about a state where synchronization
is assumed among all of the nodes within each cluster. This linearized equation is
the variational equation and it determines the stability of the clusters.
62
bc
a
0 symmetries
(11 clusters)
32 symmetries
(5 clusters)
5,760 symmetries
(3 clusters)
A =  
B =  
A =  A =  
B =  B =  
Figure 4.2: Three randomly generated networks with varying amounts
of symmetry and associated coupling matrices. (a) Nodes of the same
color are in the same synchronization cluster. The colors show the maxi-
mal symmetry the network dynamics can have given the graph structure.
(b) A graphic showing the structure of the adjacency matrices of each
network (black squares are 1, white squares are 0). (c) Block diagonal-
ization of the coupling matrices A for each network. Colors denote the
cluster, as in (a). The 2 × 2 transverse block for the 32 symmetry case
comes from one of the IRRs being present in the permutation matrices
two times. The matrices for the 32 symmetry case are given in Equations
4.13,4.14 and 4.15.
63
4.4 Master stability function formalism
The master stability function (MSF) [27] is widely used to analyze the stability
of global synchronization in networks with Laplacian coupling. My goal here is to
give a very brief review of this technique because the linear stability analysis of
cluster synchronization we present here is an extension of the MSF formalism, and
a lot of concepts from the MSF formalism will carry over to the analysis of cluster
synchrony. In particular, the concepts of the synchronization manifold and the
transverse directions will be very important in both.
We consider a system of coupled oscillators described by
x˙i(t) = F(xi(t)) + σ
∑
j
LijH(xj), i = 1, ..., N, (4.3)
This is the same as Equation 4.2, except that the Laplacian is used in place of the
adjacency matrix. By design, Equation 4.3 admits a globally-synchronous solution
xs = x1 = x2 = . . .xs. Because the coupling is Laplacian, i.e. the row sums are 0,
in the case of global synchronization we have
∑
j LijH(xs) = 0, and Equation 4.3
reduces to
x˙s(t) = F(xs(t)). (4.4)
That is, all of the oscillators evolve in time in exactly the same way as a single
oscillator would. This synchronous solution always exists, but is not always stable.
To analyze the stability of global synchrony, we consider small perturbations δxi(t)
away from the synchronous solution xs.
xi(t) = xs(t) + δxi(t) (4.5)
64
We will then expand Equation 4.3 to first order in this perturbation, obtaining an
non-autonomous, linear system of differential equations for δxi(t).
δx˙i(t) = DF(xs(t))δxi(t) + σ
∑
j
LijDH(xs)δxj(t) (4.6)
Where DF and DH are the Jacobian matrices of F and H.
An important simplification comes from expressing Equation 4.6 in the basis
of the eigenvectors of the Laplacian matrix, vi.
η i = v
T
i ⊗ Inδx (4.7)
Here, In is the n × n identity matrix, and the ⊗ symbol denotes the Kronecker
product. η i, is a vector quantity with the same dimensionality as the dynamics of a
single node, n.
η˙ i(t) = [DF(xs(t)) + σλiDH(xs)]η i(t) (4.8)
Here, i indexes the eigenvectors and eigenvalues (λi) of the Laplacian.
The strength of the MSF formalism is that it decouples network connectivity
from the dynamics. Given a particular dynamical system, (e.g. Lorenz oscillators),
the stability of global synchronization can be determined from the eigenvalues of
the Laplacian matrix alone. By writing the variational equations is this basis, the
different modes of the network dynamics uncouple. Equation 4.6 is nN -dimensional:
it has N nodes, and each node has n dimensional dynamics. In contrast, Equation
4.8 is n-dimensional. Furthermore, σλi can be viewed as a parameter that encap-
sulates the coupling strength and the network topology. The magnitudes of the
perturbations ηi(t), will either grow or shrink exponentially in time. This rate of
65
growth or decay is referred to as a transverse Lyapunov exponent. The systematic
dependence of the transverse Lyapunov exponents on the σλi is called the master
stability function.
Laplacian matrices have always a trivial eigenvector, v1 = [1, 1, . . . , 1], with
an eigenvalue of 0. In the context of synchronization, this vector has a physical
interpretation: perturbations in this direction effect all nodes equally, and so do not
break synchrony. The other eigenvectors are perpendicular to v1, and perturbations
along these directions will tend to drive the network out of synchrony.
The analysis of cluster synchronization will proceed along very similar lines. It
will not be possible to reduce the dimensionality of the problem as it is in the case of
global synchrony. In the case of cluster synchrony, there will be more synchronous
directions. If there are M clusters, there will be M orthogonal vectors that do
not break cluster synchrony. It will also be possible to find directions that are
transverse to a particular cluster. Section 4.6 describes how to find these directions,
the basis vectors of the IRR coordinate system. Section 4.5 describes the form that
the variational equations take in this coordinate system.
4.5 Variational equations in IRR Coordinates
Equation (4.2) is expressed in the “node” coordinate system, where the sub-
scripts i and j are identified with enumerated nodes of the network. Beyond iden-
tifying the symmetries and clusters, group theory also provides a powerful way to
transform the variational equations to a new coordinate system in which the trans-
66
formed coupling matrix B = TAT−1 has a block diagonal form that matches the
cluster structure, as described below. The transformation matrix T is not a simple
node re-ordering, nor is it an eigendecomposition of A. The process for computing
T is nontrivial, and involves finding the irreducible representations (IRR) of the
symmetry group. We call this new coordinate system the “IRR coordinate system.”
In section 4.6 we show the steps necessary to obtain the symmetries, clusters, and
the transformation T .
Once we have T we can transform the variational equations as follows. Let
Cm be the set of nodes in the mth cluster with synchronous motion sm(t). Then
the original variational equations about the synchronized solutions are (in vectorial
form and in the node coordinates),
δx˙(t) =
[
M∑
m=1
E(m) ⊗DF(sm(t)) + σA
M∑
m=1
E(m) ⊗DH(sm(t))
]
δx(t), (4.9)
where the Nn-dimensional vector δx(t) = [δx1(t)T , δx2(t)T , ..., δxN(t)T ]T and E(m)
is an N -dimensional diagonal matrix such that
E(m)ii =



1, if i ∈ Cm,
0, otherwise,
(4.10)
i = 1, ..., N . Note that
∑M
m=1E
(m) = IN , where IN is the N -dimensional identity
matrix.
Applying T to Eq. (4.9) we arrive at the variational matrix equation shown
in Eq. (4.11),where η(t) = T ⊗ In δx(t), J (m) is the transformed E(m), and B is the
67
block diagonalization of the coupling matrix A,
η˙(t) =
[
M∑
m=1
J (m) ⊗DF(sm(t)) + σB ⊗ In
M∑
m=1
J (m) ⊗DH(sm(t))
]
η(t), (4.11)
where we have linearized about synchronized cluster states {s1, ..., sM}, η(t) is the
vector of variations of all nodes transformed to the IRR coordinates and DF and
DH are the Jacobians of the nodes’ vector field and coupling function, respectively.
We note that this analysis holds for any node dynamics, steady state, periodic,
chaotic, etc.
We can write the block diagonal B as a direct sum
⊕L
l=1 Id(l)⊗C
l, where C l is
a (generally complex) pl× pl matrix with pl = the multiplicity of the lth IRR in the
permutation representation {Rg}, L = the number of IRRs present, and d(l) = the
dimension of the lth IRR, so that
∑L
l=1 d
(l)pl = N [110, 111]. For many transverse
blocks C l is a scalar, i.e. pl = 1. However, the trivial representation (l = 1) which is
associated with the motion in the synchronization manifold has p1 = M . Note that
the vector field F can contain a self feedback term βxi as in the experiment and
other feedbacks are possible, e.g. row sums of Aij, as long as those terms commute
with the Rg. In all these cases, the B matrix will have the same structure.
Figure 4.2c shows the coupling matrix B in the IRR coordinate system for
the three example networks. The upper-left block is an M ×M matrix that de-
scribes the dynamics within the synchronization manifold. The remaining diagonal
blocks describe motion transverse to this manifold and so are associated with loss
of synchronization. Thus, the diagonalization completely decouples the transverse
variations from the synchronization block, and partially decouples the variations
68
among the transverse directions. In this way the stability of the synchronized clus-
ters can be calculated using the separate, simpler, lower dimensional ODEs of the
transverse blocks to see if the non-synchronous transverse behavior decays to zero.
4.6 Calculating the IRR transformation
Below are the steps necessary to determine the symmetries of the network, ob-
tain the clusters, find the irreducible representations (IRRs), and the most crucial
part, calculate the transformation T from the node coordinates to the IRR coordi-
nates that will block diagonalize A, since A commutes with all symmetries of the
group [112].
Using the discrete algebra software [108,109] it is straightforward to determine
the group of symmetries of A, extract the orbits which give the nodes in each
cluster and extract the permutation matrices Rg, and use the character table of
the group and the traces of the Rg’s to determine which IRRs are present in the
node-space representation of the group. Remark: This step is discussed in any book
on representations of finite groups (e.g. Ref. [107]). After this we put each Rg into
its appropriate conjugacy class. To generate the transformaion T from the group
information and we have written code on top of the discrete algebra software which
for each IRR present constructs the projection operator P (l) [107] from the node
coordinates onto the subspace of that IRR, where l indexes the set of IRRs present.
This is done by calculating,
P (l) =
d(l)
h
∑
K
α(l)K
∑
g∈K
Rg (4.12)
69
where K is a conjugacy class, α(l)K is the character of that class for the lth IRR, d
(l)
is the dimension of the lth IRR and h is the order (size) of the group. Remark:
The trivial representation (all IRR matrices=1 and α(l) = 1) is always present and
is associated with the synchronization manifold. All other IRRs are associated with
transverse directions. Next we use singular value decomposition on P (l) to find the
basis for the projection subspace for the lth IRR. Finally, we construct T by stacking
the row basis vectors of all the IRRs which will form an N ×N matrix.
We show the results of this applied to the 32-symmetry case in Fig. 4.2:
A =


















0 1 1 1 0 1 1 1 1 1 1
1 0 1 1 1 1 0 1 1 1 1
1 1 0 1 1 1 1 1 0 1 1
1 1 1 0 1 1 1 1 1 1 1
0 1 1 1 0 1 1 1 1 1 0
1 1 1 1 1 0 1 1 1 1 1
1 0 1 1 1 1 0 1 1 1 1
1 1 1 1 1 1 1 0 1 0 1
1 1 0 1 1 1 1 1 0 1 1
1 1 1 1 1 1 1 0 1 0 0
1 1 1 1 0 1 1 1 1 0 0


















(4.13)
The nontrivial clusters are the nodes [1, 8], [2, 3, 7, 9], [4, 6], [5, 10] (the
numbering of nodes matches the row and column numbers of A). The transformation
matrix is,
70
T =
1
2



















0 0 0 0 0 0 0 0 0 0 2
0 0 0 −
√
2 0 −
√
2 0 0 0 0 0
−
√
2 0 0 0 0 0 0 −
√
2 0 0 0
0 0 0 0 −
√
2 0 0 0 0 −
√
2 0
0 −1 −1 0 0 0 −1 0 −1 0 0
−
√
2 0 0 0 0 0 0
√
2 0 0 0
0 0 0 0 −
√
2 0 0 0 0
√
2 0
0 0 0 −
√
2 0
√
2 0 0 0 0 0
0 −1 1 0 0 0 −1 0 1 0 0
0 0
√
2 0 0 0 0 0 −
√
2 0 0
0 −
√
2 0 0 0 0
√
2 0 0 0 0



















(4.14)
And the block diagonal coupling matrix is, 0.75
B =



















0 −
√
2 −
√
2 0 −2 0 0 0 0 0 0
−
√
2 1 2 2 2
√
2 0 0 0 0 0 0
−
√
2 2 1 1 2
√
2 0 0 0 0 0 0
0 2 1 1 2
√
2 0 0 0 0 0 0
−2 2
√
2 2
√
2 2
√
2 2 0 0 0 0 0 0
0 0 0 0 0 −1 −1 0 0 0 0
0 0 0 0 0 −1 −1 0 0 0 0
0 0 0 0 0 0 0 −1 0 0 0
0 0 0 0 0 0 0 0 −2 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0



















(4.15)
B has a block-diagonal structure. The 5 × 5 block describes motion in the
synchronization manifold. Each of the 5 vectors acted on by this block corresponds
to a cluster. The other blocks correspond to transverse directions. Note that the
2 × 2 block couples transverse perturbations to the red and blue clusters shown in
Figure 4.1 b. This means that the red and blue clusters are intertwined. For given
parameter values, either both or neither of these clusters synchronize.
71
4.7 Experimental system
Figure 4.3a shows the optical system used to study cluster synchronization.
Light from a 1550 nm light emitting diode (LED) passes through a polarizing beam-
splitter (PBS) and quarter wave plate (QWP), so that it is circularly polarized
when it reaches the spatial light modulator (SLM). The SLM surface imparts a pro-
grammable spatially-dependent phase shift x between the polarization components
of the reflected signal, which is then imaged, through the polarizer, onto an infrared
camera [32]. The relationship between the phase shift x applied by the SLM and
the normalized intensity I recorded by the camera is I(x) = (1− cosx) /2. The
resulting image is then fed back through a computer to control the SLM.
The SLM (Boulder Nonlinear Systems P512-1550) has an active area of 7.68 ×
7.68 mm2 and a resolution of 512 x 512 pixels. The camera (Goodrich SU320KTSW-
1.7RT/RS170) has an indium gallium arsenide focal plane array (FPA) with 320 ×
256 pixels and an area of 8 x 6.4 mm2. Using a lens assembly, the SLM was imaged
onto the camera with the magnification and orientation adjusted so that each 2 ×
2 pixel area on the SLM is projected onto a camera pixel. The camera’s frame rate
was about 8 Hz and was synchronized with the SLM’s refresh rate.
Each oscillator in equation (5) corresponds to a square patch of 16 × 16 pixels
on the SLM, which is imaged onto an 8 × 8 pixel region of the camera’s FPA. The
intensity I(x) is the average of camera pixel values in this area. The same phase
shift x is applied by each of the SLM pixels in the patch. We employ a spatial
calibration and lookup table to compensate for inhomogeneities in the SLM.
72
a b
QWP
SLM
PBS
LED
camera
lens
computer
1 mm
Figure 4.3: Experimental configuration. a) Light is reflected from the
SLM, and passes though polarization optics, so that the intensity of light
falling on the camera is modulated according the phase shift introduced
by the SLM. Coupling and feedback are implemented by a computer.
b) An image of the SLM recorded by the camera in this configuration.
Oscillators are shaded to show which cluster they belong to, and the
connectivity of the network is indicated by superimposed gray lines.
The phase shifts applied by the square regions are updated according
to equation (4.16).
73
The dynamical oscillators that form the network are realized as square patches
of pixels selected from a 32 × 32 tiling of the SLM array. Figure 4.3b shows an
experimentally measured camera frame captured for one of the 11-node networks
considered earlier in Fig. 4.2 (A full video is presented in the Supplementary Movie
1). The patches have been falsely colored to show the cluster structure, and the
links of the network are overlaid to illustrate the connectivity. The phase shift of
the ith region, xi, is updated iteratively according to:
xt+1i =
[
βI(xti) + σ
∑
j
AijI(x
t
j) + δ
]
mod 2pi (4.16)
where β is the self feedback strength, and the offset δ is introduced to suppress the
trivial solution xi = 0. Eq. (4.16) is a discrete time equivalent of Eq. (4.2). Depend-
ing on the values of β, σ and δ, Eq. (4.16) can show constant, periodic or chaotic
dynamics. There are no experimentally-imposed constraints on the adjacency ma-
trix Aij, which makes this system an ideal platform to explore synchronization in
complex networks.
4.8 Experimental results
Figure 4.4 plots the time-averaged root-mean square (RMS) synchronization
error for all four of the nontrivial clusters shown in Fig. 4.3b, as a function of the
feedback strength β, holding σ constant. We find qualitatively similar results if σ is
chosen as a bifurcation parameter with β held constant. In Fig. 4.4c-e, we plot the
observed intracluster deviations (xti − x
t) for three specific values of β indicated by
the vertical lines in Fig. 4.4a-b, showing different degrees of partial synchronization
74
Δ
x R
M
S
Δx
0
0.5
1.0
0.8
0
–0.8
2
–2
2
–2
2
–2
1.5
Δx
Δx
t t t
cl
us
te
r t
ra
ns
ve
rs
e
Ly
ap
un
ov
 e
xp
on
en
t
β
0 π
0 30 0 30 0 30
π–
2
2ππ
2
3π
2
c ed
b
a
c d ecluster 1
clusters 2 & 3
cluster 4
Figure 4.4: Experimental observation of isolated and intertwined desyn-
chronization. a) Cluster synchronization error as the self feedback, β is
varied. For all cases considered, δ = 0.525 and σ = 0.6pi. Colors indi-
cate the cluster under consideration and are consistent with Fig. 4.2. b)
MLE calculated from simulation. c-e) Synchronization error time traces
for the four clusters, showing the isolated desychronization of the ma-
genta cluster and the isolated desychronization of the intertwined blue
and red clusters.
75
that can occur, depending on the parameters.
The phase shifts imparted by each oscillator xti are recorded and used to de-
termine the degree of synchronization, as shown in Fig 3. An individual oscillator’s
deviation from synchronization at a given time as shown in Fig. 3c-e is measured by
∆xti = (x
t
i − x¯
t), where x¯t denotes a spatial average of the phases of all of the nodes
within a particular cluster at time t. To quantify the degree of synchronization
within a cluster as shown in Fig. 3a, we calculate the RMS synchronization error
∆xRMS ≡
(〈
(xti − x
t)2
〉
T
)1/2
, where 〈•〉T indicates an average over a time interval
T (here taken to be 500 iterations).
Together, Fig. 4.4a and Figs. 4.4c-e illustrate two examples of a bifurcation
commonly seen in experiment and simulation: isolated desynchronization, where one
or more clusters lose stability, while all others remain synchronized. At β = 0.72pi
(Fig. 4.4c), all four of the clusters synchronize. At β = 1.4pi (Fig. 4.4d), the
magenta cluster, which contains four nodes, has split into two smaller clusters of 2
nodes each, while the other two clusters remain synchronized.
Between β = 0.72pi and β = 1.76pi, two clusters, shown in Fig. 4.2 as red and
blue, undergo isolated desynchronization together. In Fig. 4.4a, the synchronization
error curves for these two clusters are visually indistinguishable. The synchroniza-
tion of these two clusters is intertwined: they will always either synchronize together
or not at all. While it is not obvious from a visual inspection of the network that
the red and blue clusters should form at all, their intertwined synchronization prop-
erties can be understood intuitively by examining the connectivity of the network.
Each of the two nodes in the blue cluster is coupled to exactly one node in the red
76
cluster. If the blue cluster is not synchronized, the red cluster cannot synchronize
because its two nodes are receiving different input. The group analysis treats this
automatically and yields a transverse 2× 2 block in Fig. 4.2c.
The isolated desynchronization bifurcations we observe are predicted by com-
putation of the maximum Lyapunov exponent (MLE) of the transverse blocks of
Eq. (4.11), shown in Fig. 4.4. The region of stability of each cluster is predicted by
a negative MLE. While there are four clusters in this network, there are only three
MLEs: the two intertwined clusters are described by a 2-dimensional block in the
block-diagonalized coupling matrix B. These stability calculations reveal the same
bifurcations as seen in the experiment.
4.8.1 Subgroup decomposition and isolated desynchronization
The existence of isolated desynchronizations in the network experiments raises
several questions. Since the network is connected why doesn’t the desynchronization
pull other clusters out of sync? What is the relation of ID to cluster structure and
network symmetry? Is ID a phenomenon that is common to many networks? We
provide answers to all these questions using geometric decomposition of a group that
was developed in [113,114]. This technique enables a finite group to be written as a
direct product of subgroups G = H1 × ...×Hν where ν is the number of subgroups
and all the elements in one subgroup commute with all the elements in any other
subgroup. This means that the set of nodes permuted by one subgroup is disjoint
from the set of nodes permuted by any other subgroup. Then each cluster (say, Cj) is
77
permuted only by one of the subgroups (say, Hk), but not by any others. There can
be several clusters permuted by one subgroup. This is the case of the red and blue
clusters in the 32 symmetry network in Fig.4.2, because the associated Hk cannot
have a geometric decomposition, but may have a more structured decomposition
such as a wreath product [115].
We can show that the above decomposition guarantees that the nodes as-
sociated with different subgroups all receive the same total input from the other
subgroups’ nodes. Hence, nodes of each cluster will not see the effects of individ-
ual behavior of the other clusters associated with different subgroups. This enables
the clusters to have the same synchronized dynamics even when another cluster
desynchronizes. If that state is stable we have ID.
To see this let Hk, a subgroup of G, permute only cluster Cm and pi be the
permutation on the indices of nodes in Cm for one permutation Rg, g ∈ Hk. Assume
xi is not in Cm so it is not permuted by Rg and recall that A commutes with all
permutations in G, then we have (just concentrating on the terms from Cm),
[Rgx˙(t)]i = x˙i(t) = ...+ σ[RgAH(x)]i = ...+ σ[ARgH(x)]i = ...+ σ
∑
j∈Cm
AijH(xpi(j)),
(4.17)
where pi(l) is, in general, another node in Cm and the sums over other clusters are
unchanged. This shows that all nodes in Cm are coupled into the ith node in the
same way (the same Aij factor). Similarly, if we use a permution Rg′ on the cluster
Cm′ containing xi we can show that all the nodes of Cm′ are coupled in the same
78
way to the nodes in Cm. Hence, nodes of Cm′ each receive the same input sum from
the nodes of Cm whether the nodes of Cm are synchronized or not. This explains
how the cluster Cm can become desynchronized, but the nodes of Cm′ can still be
synchronized – they all have the same input despite the Cm desynchronization, thus
making the Cm′ synchronous state flow invariant. If it is also stable, this is the case
of ID. This argument is easily generalized to the case when Hk permutes nodes of
several clusters as this will just add other similar sums to Eq. (4.17). The latter case
explains the intertwined desynchronization in the experiment and is a more general
form of ID.
4.9 Generalizations
The symmetry-based analysis described allows for the linear stability analysis
of cluster synchrony in arbitrary networks. This is a major step forward, but there
is more to do.
So far, we have described how to analyze the stability of one synchronization
pattern. The orbits of the isomorphism group specify the maximally-symmetric
pattern, the pattern with the smallest number of clusters possible with adjacency
matrix coupling. We have already seen that other cluster-synchronized states are
possible. After an isolated desynchronization, the network is often in a dynamical
state that still has clusters. We need a formalism that will predict the stability of
these patterns.
We also have not considered the case of Laplacian coupling. We already know
79
that Laplacian coupling allows for global synchronization, and we will see that it
allows for more complex synchronization patterns as well.
We will see that a given network can support a large number of synchroniza-
tion patterns. One very simple example of a 5-node network we will consider has
12 possible synchronization patterns. I don’t know how the number of patterns
will scale with the size of the network, but the number of possible synchronization
patterns in a large network has the possibility of being combinatorially huge. Con-
sider, for example, how many ways an all-to-all coupled network of N oscillators
could possibly split into two equal clusters. Of course, a huge variety of unequal
clusterings are also possible in an all-to-all coupled network [116]. As I will describe
in Chapter 5, a nonlocal ring coupling configuration also supports a vast (and yet
unquantified) number of multistable “incoherent” states.
So, for now, enumerating all of the synchronization patterns possible in a large
network is a challenge, because there are often a lot of them. But, thanks to recent
work by Francesco Sorrentino and Lou Pecora, it is now possible to predict and
analyze the stability of every possible cluster-synchronous state [117]. In this section,
I will give a simple example network with 5 nodes, and 12 possible synchronization
patterns. I will describe how to predict those patterns, and briefly indicate how the
linear stability analysis will proceed. Details can be found in [117]. Finally, I will
show the experimental confirmation of some of the new types of patterns.
80
4.9.1 Subgroup clusters
The network we will consider is given by the Laplacian matrix
L =






−3 1 0 1 1
1 −3 1 0 1
0 1 −3 1 1
1 0 1 −3 1
1 1 1 1 −4






. (4.18)
This network and all of the cluster synchronization patterns it supports are shown
in Figure 4.5. The patterns A1,..., A5 arise due to subgroups of the graph’s auto-
morphism group. The decomposition of a group into a direct product of subgroups,
i.e. G = H1 × ...Hi...×Hν , can be accomplished using computational group theory
software [113, 114]. An individual subgroup Hi can be excluded, and a new group
can be constructed G ′ = H1 × ... × Hν . The new group G ′ has its own orbits, and
the orbits of each subgroup will yield a cluster pattern, just as the orbits of the
full group do. The construction of the variational equations follows the same lines
as described in Sections 4.6 and 4.5, except one considers G ′ rather than the whole
group.
4.9.2 Laplacian coupling and merged clusters
Laplacian coupling allows for merged clusters. For example, L1 is the merger of
the trivial cluster, node (5), with the cluster (1,2,3,5). In general, merged clusters are
formed by merging clusters that arise due the symmetries of the adjacency matrix.
Not all clusters can merge. For example, in the pattern A2, node (5) can’t join either
of the other clusters. There isn’t really a good way to determine which mergers
are possible, except by evaluating each one individually. It is straightforward, but
81
1 2
4
5
3
x
1
≠ x
2
≠ x
3
≠ x
4
≠ x
5
A5
1 2
4
5
3
x
1
= x
3
1 2
4
5
3
x
2
= x
4
1 2
4
5
3
x
1
= x
3
= x
5
1 2
4
5
3
x
2
= x
4
= x
5
A4
1 2
4
5
3
1 2
4
5
3
x
1
= x
3
= x
5
, x
2
= x
4
1 2
4
5
3
x
2
= x
4
= x
5
, x
1
= x
3
A3
1 2
4
5
3
x
1
= x
3
, x
2
= x
4
1 2
4
5
3
x
1
= x
4
, x
2
= x
3
x
1
= x
2
, x
3
= x
4
A2
1 2
4
5
3
x
1
= x
2
= x
3
= x
4
1 2
4
5
3
x
1
= x
2
= x
3
= x
4
= x
5
A1
L4
L3
L1
Figure 4.5: All of the synchronization patterns of a small network. The
patterns A1,...,A5 are due to symmetries of the adjacency matrix. The
patterns L1, L3, and L4 are only possible with Laplacian coupling.
82
slightly tedious, to verify by hand that the synchronization patterns L1, L3, and L4
satisfy the equations of motion.
There is a better way to determine if it is possible for a given merger to
happen than direct examination of the equations of motion. Note that, according
to Equation 4.3, the coupling term vanishes between synchronized nodes. Deleting
links between nodes that are synchronized in a proposed merged-cluster state leads
to a new, dynamically-equivalent network. If the proposed synchronization pattern
is allowed by the symmetries of the dynamically-equivalent network, then the merged
cluster state is possible. Merged cluster states will have both new synchronous and
new transverse directions.
4.9.3 Experimental observation
We show subgroup and Laplacian clusters using the experimental system de-
scribed in 4.7. Figure 4.6 shows the experimental synchronization error for each
synchronization pattern as a function of the parameter σ. For these measurements,
β = 1.45pi, and σ was decreased from pi to 0 in steps of pi/50, with 1000 iterations
at each value of σ. The phases xi were not reinitialized when the values of σ were
updated. We plot the pattern synchronization error. This quantity measures the
extent to which the nodes synchronize in one of the patterns shown in Figure 4.5,
and is defined as the average over all clusters under consideration (e.g. the clusters
(1, 5, 3) and (2, 4) for L3), of the time average of 〈|xi − 〈x〉 |〉, where 〈•〉 denotes the
average over the nodes in a cluster.
83
0
π π 3π
π
424
0
1
2
3
4
5
S
y
n
c
h
ro
n
iz
a
ti
o
n
 e
rr
o
r
S
ta
b
il
it
y
A4
L4
A3
L3
A2
A1
L1
σ
Figure 4.6: The experimental synchronization error for each synchro-
nization pattern as a function of the parameter σ. Underneath we plot
the results of our stability analysis applied to each one of the cluster
synchronization patterns, where a colored dot labels the values of σ for
which the corresponding pattern is stable.
84
Chapter 5: Chimeras in Coupled-Map Lattices
This chapter is based on work from the following publication
Aaron M. Hagerstrom, Thomas E. Murphy, Rajarshi Roy, Philipp Ho¨vel, Iryna
Omelchenko, and Eckehard Scho¨ll. Experimental observation of chimeras in coupled-
map lattices. Nature Physics, 8(9):658–661, 2012.
5.1 Introduction and outline
Chimera states were first recognized in a ring of coupled phase oscillators by
Kuramoto and Battogtokh in 2002 [29]. Networks of phase oscillators are widely
studied in part because the phase oscillator model describes the generic behavior
of systems near a Hopf bifurcation with weak coupling [28], and in part because
in many cases low-dimensional models of infinitely large populations can be con-
structed [118]. Depending on the strength of the coupling, a large population of
globally-coupled phase oscillators will either show synchronous “coherent” behavior,
or unsynchronized “incoherent” behavior. When nonlocal coupling is used instead,
both behaviors can be present at the same time: networks with nonlocal coupling
can have large, coexisting spatial regions of synchrony and incoherence. These re-
gions of qualitatively different behaviors are surprising because the coupling in these
models is homogeneous. There’s no reason why one particular part of a ring should
85
synchronize while another part does not. This coexistence of seemingly incongru-
ous behaviors inspired Abrams and Strogatz to name these dynamical states after
a creature from Greek mythology with a goat’s head, a lion’s head, and a snake for
a tail [30]. Hundreds of papers have been written about Chimera states exploring
various types of coupling, different local dynamics, coupling with time delay, the ef-
fect of parameter heterogeneity, their metastability as a function of size, and many
other things [30,119–121,121–129].
Before 2012, there were no experimental demonstrations of chimeras. Now
there are several. This chapter describes one of the first two experimental real-
izations of chimera states. The other experiment was performed by Mark Tinsley,
Simba Nkomo, and Ken Showalter, who observed chimera states in two coupled pop-
ulations of chemical oscillators [33]. Their experiment was inspired by a theoretical
paper from 2008 which presented a low-dimensional analytical model of two coupled
populations of phase oscillators [119]. This theoretical paper was also the inspiration
for an experiment by Erik Martens, Shashi Thutupalli, Antoine Fourrie`re, and Oskar
Hallatschek, who used metronomes as phase oscillators, and constructed two coupled
populations by placing the metronomes on swings coupled by a spring [36]. Chimera
states have also been observed in time-delayed feedback loops [37,38]. Time-delayed
systems have an infinite-dimensional phase space which can be thought of as analo-
gous to a spatial domain. Chimera states have also been observed in electrochemical
systems including some with global coupling [34,35].
Many models which display chimera states have nonlocal coupling. In a spatial
domain, oscillators are said to be nonlocally coupled if each oscillator has a finite
86
range of influence, rather than being coupled to the whole system, or only to an
infinitesimally small local area. In systems consisting of a finite number of oscillators,
coupling is said to be nonlocal if a given oscillator is coupled to other oscillators
within a distance that is an appreciable fraction of the system size. Global coupling,
in which all oscillators receive the same signal [116], and local coupling [130–135]
were much more widely studied prior to the discovery of chimeras, in part because it
was hard to think of experimental situations in which nonlocal coupling might come
about. Shima and Kuramoto argued that nonlocal coupling could arise in reaction-
diffusion systems if one reagent diffused much more slowly than another [31].
When we wrote our paper in 2012, it was generally believed that chimera states
require nonlocal coupling. Now, in 2015, chimeras have been described in models
with purely local coupling [136] and in both models [137] and experiments [34, 35]
with global coupling. It’s also worth noting that more than a decade before chimeras
in phase oscillators were described by Kuramoto and Battogtokh, Kaneko described
dynamical states in globally-coupled iterated maps with a coexistence of a few large
synchronized clusters and many small unsynchronized clusters [116]. These states
probably deserve to be called chimeras.
The chimera states we describe here exist in a coupled-map lattice (CLM)
[130–135] in one or two [138] dimensions. Although we use the term “oscillator” to
describe the iterated maps in the experiment, the local dynamics of an individual
oscillator are not oscillatory. Each oscillator acts like a chaotic iterated map when
it is not coupled to other oscillators. Even though individual oscillators behave
chaotically, we will see situations where a network of coupled oscillators is periodic.
87
The network is updated in discrete time, and so the temporal periodicity of a given
lattice site will be an integer number of iterations. In the chimeras we describe,
all of the nodes are synchronized, in the sense that the entire system is periodic
with the same period. In this case the terms “coherence” and “incoherence” refer
to spatial behavior rather than temporal behavior. We use the term coherent to
refer to smooth profiles, and incoherent to refer to spatial profiles with a large
number of discontinuities. These types of chimeras in coupled map lattice were first
theoretically considered and numerically simulated in refs. [139,140].
5.1.1 Organization of this chapter
Section 5.2 describes the chimera states we observed, and the dependence of the
dynamics on parameter values in both simulation and experiment. Two parameters
are particularly important: the coupling radius r, and the coupling strength .
At large , we see smooth spatial profiles. At low , we see incoherent patterns
with many spatial discontinuities. Section 5.3 derives an approximation for a critical
coupling strength at which a discontinuity develops.
Our chimera states have incoherent regions. These should be thought of as
finite-sized regions, rather than thin boundaries between neighboring smooth re-
gions. In Section 5.4, we show simulation data which demonstrates that the inco-
herent region scales with the number of oscillators, N .
We see patterns with different spatial wavenumbers. Intuitively, a smaller
coupling radius leads to patterns with higher wavenumber. In Section 5.5 we show
88
that patterns with different wavenumbers are related by a simple scaling argument.
5.2 Coupled-map lattices in experiment
Figure 5.1 shows the experimental setup of the optical CML. This is the same
device described in Chapter 4. The device employs a phase spatial light modulator
(SLM) and an infrared camera in feedback configuration. Both the SLM and the
camera frames are partitioned into an M×M array of square regions. These regions
correspond to nodes in the network of coupled maps. The phase shift xi,j and
intensity I for a given region (i, j) are related by
I(xi,j) =
1
2
(1− cosxi,j) , i, j = 0, . . . ,M − 1. (5.1)
The intensity has been normalized to be between 0 and 1. A lens is used
to project an image of the SLM onto an infrared camera, which records the two-
dimensional intensity pattern I(xi,j). We construct a network of iterated maps
by using the computer to communicate between the camera and the SLM. Time
evolution of the network is achieved by iteratively updating the phase applied by
each region of the SLM in a way that depends upon the intensity measured by the
camera.
5.2.1 Chimeras in 1D and 2D
We present results for two different coupling schemes shown schematically in
Figures 5.1(b) and 5.1(c). In the one-dimensional (1D) configuration, the elements
in the array are arranged as a ring with periodic boundary conditions. The SLM is
89
spatial light
quarter-wave plate
polarizing
beamsplitter
1550 nm LED
imaging optics
coupling and feedback
      (through a computer)
camera
modulator (SLM)
… 
…
……
0,M-1 
M-1,M-1
 
M-1,0 
0,0 0,1 0,2 0,3 0,4 
1,0 1,1 1,2 1,3 1,4 
2,0 2,1 2,2 2,3 2,4 
3,0 3,1 3,2 3,3 3,4 
4,0 4,1 4,2 4,3 4,4 
 one-dimensional configuration
0 1 2 3 4 5 6 N-1 N-2 
R 
… 
two-dimensional configuration
a)
b) c)
        
 
  
Figure 5.1: Experimental device and coupling configuration. a) Op-
tical configuration. Polarization optics create a nonlinear relationship
between the spatially dependent phase shift applied by the SLM and
the intensity of the light falling on the camera. Feedback and coupling
are implemented using a computer. b) One- and c) two-dimensional
coupling configurations are shown schematically. Because the elements
are coupled diffusively to their neighbors within a range R in either one
or two dimensions with periodic boundary conditions, the coupling is
identical for all oscillators.
90
treated as a 1D lattice with the elements coupled in a raster-ordered arrangement.
If xti is the phase of the ith element in the ring at the tth iteration, and I(x
t
i) is
the intensity measured by the corresponding region in the camera, then the phase
is updated according to
xt+1i = 2pia
{
I
(
xti
)
+

2R
R∑
j=−R
[
I
(
xti+j
)
− I
(
xti
)]
}
(5.2)
where i = 0, . . . , N−1 in the 1D ring configuration and the index i is periodic modulo
N . Thus each element is coupled diffusively to all of the elements within a distance
R on either side, and  describes the strength of the coupling. The parameter
a controls the temporal dynamics of an isolated map. We choose a = 0.85 such
that the local map xt+1 = f(xt) ≡ 2piaI(xt) is chaotic with a Lyapunov exponent
λ ≈ 0.58.
We also examine a two-dimensional (2D) coupling scheme. In this configu-
ration, each element is coupled to its neighbors within a square region, and the
boundary conditions are periodic in both dimensions.
xt+1i,j = 2pia
{
I
(
xti,j
)
+

4R2
R∑
k,l=−R
[
I
(
xti+k,j+l
)
− I
(
xti,j
)]
}
(5.3)
Figures 5.2(a) and 5.2(b) characterize the dependence of the dynamics of the
1D and the 2D systems on coupling strength and range by plotting a measure of
temporal entropy (see Section 5.2.2) in order to characterize the periodicity of the
dynamics. The color maps were obtained through numerical simulation of equations
(5.2) and (5.3). A period-2 profile will have an entropy of 1 bit, and a chaotic profile
will have an entropy close to 9 bits, since 512 bins were used. We see spatial profiles
91
a)
c)
b)
d)
0 1
S
im
ul
at
io
n
E
xp
er
im
en
t
Coupling radius r
C
ou
pl
in
g 
st
re
ng
th
 ε En
tr
op
y
9bits
1 bit
i/N
0
1
0.0 0.1 0.2 0.3 0.4
0.0
0.2
0.4
0.6
0.8
1.0
0.5
j/N
In
te
ns
ity
1.00
0.25
D
D
E
E F
F
D
E
F
K=1
Incoherence
Chaos
Global Sync.
0 1
S
im
ul
at
io
n
E
xp
er
im
en
t
Coupling radius r
C
ou
pl
in
g 
st
re
ng
th
 ε
In
te
ns
ity
i/N
0.25
1.00
0.0 0.1 0.2 0.3 0.4 0.5
0.0
0.2
0.4
0.6
0.8
1.0
B
B
A
A
C
C
A
B
C
K=1
Incoherence
Chaos
Global Sync.
Figure 5.2: Parameter space of the 1D and 2D CML systems. a) Pa-
rameter space of the one-dimensional configuration (N = 256 elements).
Blue corresponds to low entropy and periodicity, orange corresponds to
high entropy and chaos. The dashed line indicates the critical coupling
strength c = 0.54 b) Parameter space of the two-dimensional config-
uration (128 × 128 lattice). There is a prominent K = 1 tongue, and
two tongues with more complex spatial patterns. c) Experimental and
numerical realizations of the 1D system. In the middle panel, B, the inco-
herent region is highlighted in yellow. Labels A (coherence), B (chimera),
C (incoherence) show positions in the parameter space of panel a). d)
Experimental and numerical realizations of the 2D system. In the middle
panel, E, the incoherent region is highlighted.
92
that are periodic in time, as well as profiles that are chaotic in time. There is a rich
variety of both coherent and incoherent spatial structures. The thin vertical regions
marked K = 1, 2, 3 in figure 5.2(a) indicate conditions for which the dynamical
state is periodic in time and forms a spatially periodic standing-wave pattern with
a spatial wavenumber of K. Within each of these regions, there is a transition from
spatial coherence to incoherence with decreasing coupling strength , as predicted
in [139,140].
Equations (5.2) and (5.3) admit a globally synchronous solution. The hatched
region in Figures 5.2(a) and 5.2(b) indicate the region in parameter space for which
global synchrony is stable, which is identical between the 1D and 2D systems.
A comparison of the 1D and 2D systems reveals that these systems show equiv-
alent behavior within a limited region in parameter space. In the 2D system there
is a prominent region marked K = 1, which contains profiles that are homogeneous
in one direction. Every solution of the 1D problem corresponds to a solution of the
2D problem that is uniform in one direction, i.e., xi,j = xi. While K = 2, 3 . . .
profiles also exist in principle, only the K = 1 profile is observed in simulation and
experiment in the 2D system within the evolution time we consider. Simulations
also reveal two smaller regions that contain complex profiles with spatial structures
that cannot be realized in one dimension.
Figure 5.2(c) shows experimental and numerical realizations of three profiles
from the K = 1 region. The sequence A,B,C shows a progression from spatial
coherence to spatial incoherence. All profiles have a temporal period of 2. At a
coupling radius r ≡ R/N = 70/256 ≈ 0.27 and  = 0.75 (A), the profile has a
93
smooth spatial variation. At r = 105/256 ≈ 0.41,  = 0.44 (B), the profile shows
two large domains of coherent, synchronized behavior separated by narrow but finite,
regions of incoherent behavior highlighted in yellow. Numerical evidence indicates
that the width of these regions as a fraction of the system size remains constant as N
increases, and hence this is not a boundary effect (see Section 5.4). While the entire
profile is periodic in time, the dynamics of the coherent and incoherent regions are
qualitatively different in their degree of spatial coherence: within the incoherence
intervals there are numerous admissible combinations of upper and lower states
whose multiplicity scales exponentially with the system size [139]. This mixture
of qualitatively different behaviors is analogous to the chimera states discussed in
references [29,30,119,122–124], where an array of identical oscillators splits into two
domains: one coherent and phase locked, the other incoherent and desynchronized.
Unlike the chimeras in continuous-time phase oscillators [29, 30], the generalized
chimeras here and in [139,140] have two coherent parts (with high and low intensity,
respectively) and two incoherent parts, as a result of their mechanism of nascence
from the completely coherent K = 1 spatial profile. Finally, at r = 115/256 ≈ 0.45,
 = 0.375 (C), the profile is completely incoherent. The same scenario occurs in
the K = 2 and K = 3 regions and can be interpreted through a spatial rescaling
described in Section 5.5.
Figure 5.2(d) shows snapshots of experimental and numerical realizations of
the 2D system in a 32× 32 lattice, obtained with similar parameter values to those
used in Figure 5.2(b). All of these realizations are periodic in time with period 2,
except (E), which has a period of 4 in the experimental realization shown. Like
94
the 1D system, the 2D system undergoes a coherence-incoherence transition as  is
decreased. At r = 9/32 ≈ 0.28,  = 0.75, the system exhibits a smooth profile (D).
For r = 13/32 ≈ 0.41 and  = 0.44, there is a chimera-like coexistence of coherence
and incoherence (E). Finally, at r = 14/32 ≈ 0.44 and  = 0.375, the system is fully
incoherent (F).
5.2.2 Periodicity quantification
In Figures 5.2(a) and 5.2(b), the color of each point corresponds to a single
numerical simulation. Lattice sizes of 256 for the 1D case and 128× 128 for the 2D
case were used. Initially, the phase of each node is random and uniformly distributed
between 0 and 2pi. In each simulation we discard 50000 transient iterations. Using
the next 5000 iterations, we construct one histogram for each of the N nodes in the
network, binning the values the phase achieves in these iterations. 512 bins were
used. Thus we estimate pni , the fraction of the time that node i will spends in the
nth bin. From these statistics we obtain the entropy, which is then averaged over
the N nodes in the system
H = −
1
N
N−1∑
i=0
511∑
n=0
pni log2 p
n
i . (5.4)
We note that for this system the temporal behavior of the whole network is either
chaotic or periodic, and the spatial average of the entropy will not differ substan-
tially from the entropy of any given node. We also note that this calculation only
characterizes the temporal behavior of the system, and does not distinguish between
spatial coherence and incoherence.
95
5.3 Critical coupling strength
The development of a spatial discontinuity can be explained analytically. We
derive the critical coupling strength c at which the smooth profile breaks up giving
rise to the incoherence regions (B) of figure 5.2(c). Considering the spatially contin-
uous version of equation (5.2) and solutions with wavenumber K = 1 and period-2
dynamics in time, the system dynamics are given by alternating profiles x0(u) and
x1(u) for even and odd time steps, respectively. By evaluating the spatial derivative
of this equation with the profiles x0(u) at the positions where the smooth profiles
break up, we obtain two equations which, when multiplied, will yield a condition
for the critical coupling strength.
Starting from the system equation (5.2)
xt+1i = f
(
xti
)
+

2R
R∑
j=−R
[
f
(
xtj+i
)
− f
(
xti
)]
, (5.5)
where we define the local map as f(x) = 2piaI(x), we obtain its spatially continuous
version for N →∞ and r = R/N :
xt+1(u) =f
(
xt(u)
)
+

2r
∫ u+r
u−r
[
f
(
xt(v)
)
− f
(
xt(u)
)]
dv. (5.6)
Let us consider a solution with wave number K = 1 and period-2 dynamics in time.
Hence we can reduce the dynamics by even and odd time steps x0(u) and x1(u),
respectively:
x1−j(u) = (1− )f
(
xj(u)
)
+

2r
u+r∫
u−r
f
(
xj(v)
)
dv (5.7)
96
-2
-1
 0
G
a)
 3
 5
0 255Index i
b) x
x
i
i
0
1
u1 u2
x∗
 0.25
 1
0 N-1
Inte
nsi
ty
Index i
e) 0.25
 1
Inte
nsi
ty
i’=0 i’=N-1 i’=2N-1d)
 0.25
 10 3N-1
Inte
nsi
ty
Rescaled index i’c)
Figure 5.3: Critical coupling and scaling of coherent profiles. a) Condi-
tion G for breaking of the smooth profile. The solid curve corresponds to
the critical coupling strength c = 0.54. The dashed and dotted curves
refer to  = 0.675 > c and  = 0.475 < c, respectively. b) Snap-
shots of the period-2 solution for even (x0i (black)) and odd (x
1
i (red))
time steps at the critical coupling strength of panel (a). The light red
and gray dots refer to the experiment. Points u1, u2 are discontinuity
points, where the profile breaks, and x∗ marks the fixed point of the local
map. Other parameters: N = 256, a = 0.85, r = 88/256 ≈ 0.34. c)-e)
Scaling of coherent intensity profiles. The measured coherent profiles
with spatial periods K = 1, 2, 3 shown in black for coupling radius (c)
r3 = 22/256 ≈ 0.09, (d) r2 = 33/256 ≈ 0.13, and (e) r1 = 66/256 ≈ 0.26.
The profile K = 1 numerically obtained from equations (5.2),(5.1) and
its rescaled profile are shown in red in panels c) and d),e), respectively.
Other parameters: a = 0.85,  = 0.8, and N = 256.
97
with j = 0, 1. Taking the spatial derivative yields
∂ux
1−j(u) =(1− )∂uf
(
xj(u)
)
∂ux
j(u)+
+

2r
[
f
(
xj(u+ r)
)
− f
(
xj(u− r)
)]
. (5.8)
At the points where the smooth profile breaks up, the spatial derivative ∂uxj(u)
becomes infinite. Considering that ∂ux0(u), ∂ux1(u) diverge to infinity, we can ne-
glect the coupling term on the right-hand side. Multiplying the equations for even
and odd time steps we obtain:
∂ux
0(u)∂ux
1(u) =
[
(1− )2∂uf
(
x0(u)
)
∂uf
(
x1(u)
)]
× ∂ux
0(u)∂ux
1(u), (5.9)
which yields the following condition
1 = (1− )2∂uf
(
x0(u)
)
∂uf
(
x1(u)
)
. (5.10)
As the local dynamics is governed by the map f(x) = pia(1−cosx) and its derivative
is equal to f ′(x) = pia sinx, we introduce the function
G(u) = (pia)2(1− )2 sin(x0(u)) sin(x1(u))− 1 (5.11)
as a deviation from the previous condition (5.10).
For fixed coupling radius r and fixed  we can numerically calculate x0(u), x1(u).
For the critical coupling strength c the condition G(u∗) = 0 holds at discontinuity
points u∗.
In order to obtain an analytic approximation for the critical coupling strength
c, we assume x0(u) = x1(u) = x∗ at the discontinuity points u∗, where x∗ is a fixed
98
 0
2π
2π
x
f
f2
x*
0
x*
1
x*
2
Figure 5.4: Local map f(x) = 2piaI(x) = pia(1 − cosx) (red) and its
twice iterate f 2 (blue) for a = 0.85. The filled and open circles refer to
the fixed points of f and the twice iterated f 2.
99
point x∗ = pia(1− cosx∗) of the map, hence
G(u) = (pia)2 (1− c)
2 sin2 x∗ − 1 = 0. (5.12)
This leads to
c = 1∓
1
(pia)| sinx∗|
, (5.13)
where the minus sign should be chosen since the lower value of c represents the
threshold where the smooth profile breaks up with decreasing coupling strength ,
and the spatial coherence is lost. The map f(x) = pia(1 − cosx) has three fixed
points x∗0 = 0, x
∗
1 ≈ 0.79, and x
∗
2 ≈ 4.13 for a = 0.85, as shown in Figure 5.4.
The solutions with K = 1 in space and period-2 in time shown in Figure 5.3(b)
can be understood qualitatively as follows: part of the network is in the upper state,
part is in the lower state, and both states alternate with period 2. These two states
correspond very roughly to the two fixed points of the twice iterated map above and
below the fixed point x∗2 ≈ 4.13. They have bifurcated from x
∗
2 which has become
unstable in the twice iterated map, and the discontinuity point between these is the
fixed point x∗2 (see Figure 5.4). Hence, the approximation x
∗ ≈ x∗2 in Equation (5.13)
leads to c ≈ 0.55, which is close to the result obtained numerically in Figures 5.3(a)
and 5.3(b).
One can even obtain an approximate full analytical solution (for a > 0.6) by
Taylor expanding the local map f about pi using
cosx∗2 ≈ −1 +
1
2
ψ2, (5.14)
100
where ψ = x∗2 − pi. This gives an equation for the fixed point
pi + ψ ≈ pia
(
2−
1
2
ψ2
)
(5.15)
and hence
ψ ≈
1
pia
(−1 +
√
1 + 2pi2a(2a− 1)). (5.16)
This gives ψ ≈ 0.96, x∗2 ≈ 4.10, and c ≈ 0.54, which agrees with the observed
behavior of both the experimental system and simulations to a few percent.
5.4 Scaling of the incoherent regions in chimera states with N
To illustrate the relationship between the size of the incoherent regions high-
lighted in yellow in Figs. 5.2(c),(d), and the finite size (N) of the experimental
lattice, we focus on the 1D case in Fig. 5.2(c) panel (B). While these incoherent
regions are small, they are finite and not just a boundary between the coherent
domains. The size of the incoherent region as a fraction of the system size remains
approximately constant as the system size is increased in numerical simulations,
which indicates that it is not a finite-size boundary effect. This is shown in Fig. 5.5,
which depicts snapshots for N = 256 and N = 1024, and a decreasing sequence of
coupling strengths  = 0.44, 0.42, and 0.40, with the same coupling radius r ≈ 0.41
as in panel B of Fig. 5.2(c). Furthermore, Fig. 5.5 demonstrates that the inco-
herence regions become wider for decreasing coupling strength  as expected in
the transition scenario from spatial coherence to complete incoherence via chimera
states (cf. [139, 140]). The initial conditions in Fig. 5.5 are carefully prepared to
101
 0.4
 1
 0  255Int
e
n
s
i
t
y
Index i
(a) ε=0.44
 0.4
 1
 0  255Int
e
n
s
i
t
y
Index i
(c) ε=0.42
 0.4
 1
 0  255Int
e
n
s
i
t
y
Index i
(e) ε=0.40
 0.4
 1
 0  1023Int
e
n
s
i
t
y
Index i
(b) ε=0.44
 0.4
 1
 0  1023Int
e
n
s
i
t
y
Index i
(d) ε=0.42
 0.4
 1
 0  1023Int
e
n
s
i
t
y
Index i
(f) ε=0.40
 0.4 1
 50  75Intensit
y
Index i  0.4 1
 40  90Intensit
y
Index i  0.4 1
 20  100Intensit
y
Index i
 0.4 1
 200  300Intensit
y
Index i  0.4 1
 160  360Intensit
y
Index i  0.4 1
 80  400Intensit
y
Index i
Figure 5.5: Simulations to illustrate scaling with system size of the in-
coherent regions for the 1D system. Snapshots for (a),(c),(e) N = 256,
and (b),(d),(f) N = 1024. Decreasing coupling strength is used in panels
(a),(b)  = 0.44, (c),(d)  = 0.42, (e),(f)  = 0.40. Other parameters:
a = 0.85, r = 0.41. Transients of 1000 time steps were neglected. The
insets show the enlarged left incoherence region for each snapshot.
generate chimera states. We divided the ring of N elements into four equal parts,
and used I(x0i = 4.5) ≈ 0.6 and I(x
0
i = 2.5) ≈ 0.9 on the segments
i
N
∈
[
−
1
8
,
1
8
]
and
i
N
∈
[
3
8
,
5
8
]
, respectively, and random sequences of these two values for the
other two parts, which contain the incoherent domains as subsets [139, 140]. The
experimental measurements presented in Fig. 5.2 were made with random initial
conditions. The role of initial conditions in the emergence of chimera states and the
existence of multistability needs to be explored in the future.
102
 3
 5
0
x i
i
a)
 3
 5
0 N-1
0 2N-1
i
i’b)
 3
 5
0 N-1
0 3N-1
i
i’c)
x i
x i
Figure 5.6: Coherent profiles with space periods K = 1, 2, 3 accord-
ing to equation (5.2) for coupling radius (a) r1 = 66/256 ≈ 0.26, (b)
r2 = 33/256 ≈ 0.13, and (c) r3 = 22/256 ≈ 0.09 shown in black. The
rescaled coherent profile K = 1 in panels (b),(c) is shown in red. Other
parameters: a = 0.85,  = 0.8, and N = 256.
103
5.5 Spatial rescaling and coherent patterns with K=1,2,3...
To obtain the scaling relation for coherent profiles, we consider again the
thermodynamic limit N →∞ and the spatially continuous version of system (5.2):
xt+1(u) = f(xt(u)) +

2r
u+r∫
u−r
[
f(xt(v))− f(xt(u))
]
dv, (5.17)
where f(x) = 2piaI(x), r = R/N is a coupling radius, and the coherent solution
xti, i = 0, . . . , N − 1, of the system (5.2) approaches a corresponding smooth profile
xt(u), where u ∈ [0, 1].
In the following, we consider the case that the coupling strength  is fixed.
Let there exist a state xK′(u) with spatial period K ′. The dynamics of this state is
governed by
xt+1K′ (u) = f(x
t
K′(u)) +

2r
u+r∫
u−r
[
f(xtK′(v))− f(x
t
K′(u))
]
dv, (5.18)
The solution xK′(u) will be compared with the state xK(u) with spatial period K.
Our goal is to show that the dynamics of xK′(u) is identical to the dynamics of xK(u˜)
at a rescaled space position u˜ = (uK ′)/K for a correspondingly rescaled coupling
radius r˜ = (rK ′)/K. For simplicity, we will consider that K ′ = 1, and we assume
that xK′(u) = xK(u/K). The dynamics of the solution xK of Equation (5.18) with
rescaled u→ u/K and rescaled coupling radius r → r/K is governed by
xt+1K
( u
K
)
= f
(
xtK
( u
K
))
+
K
2r
(u+r)/K∫
(u−r)/K
[
f
(
xtK(v)
)
− f
(
xtK
( u
K
))]
dv, (5.19)
Let v = w/K, then
xt+1K
( u
K
)
= f
(
xtK
( u
K
))
+

2r
u+r∫
u−r
[
f
(
xtK
(w
K
))
− f
(
xtK
( u
K
))]
dw. (5.20)
104
With the use of the above mentioned assumption xtK(u/K) = x
t
K′(u), we obtain
xt+1K
( u
K
)
= f(xtK′(u)) +

2r
u+r∫
u−r
[
f(xtK′(w))− f(x
t
K′(u))
]
dw = xt+1K′ (u), (5.21)
i.e., identical dynamical evolution.
This is depicted in Figure 5.6, which illustrates the scaling of the coherent
profiles obtained from numerical simulations of equation (5.2) with a = 0.85,  =
0.8, and N = 256. Coherent profiles with K = 1, 2, 3 for a coupling radius (a)
r1 = 66/256 ≈ 0.26, (b) r2 = 33/256 ≈ 0.13, and (c) r3 = 22/256 ≈ 0.09 are shown
in black. The rescaled coherent profile for K = 1 in the panels (b) and (c) is shown
in red.
To conclude this argument, if in the system (5.17) the coupling strength 
is fixed and there exist two solutions xK′ and xK with space periods K ′ and K,
respectively, then xK′(u) = xK (uK ′/K) for corresponding coupling radii r and
(rK ′)/K, respectively.
105
5.6 Concluding remarks
In this chapter, we showed that the CLM chimeras described in [139,140] can
be implemented in experiment. We also described some of their basic properties,
including the development of a spatial discontinuity with decreasing , and the
existence of patterns with different spatial wavenumbers. We also showed that
there are both one-dimensional and two-dimensional systems with similar behavior,
although, strangely, Chimeras with wavenumbers K other than 1 are not seen in
the two-dimensional system. The multistability of the chimera states, however, is
an important question that remains unaddressed.
As far as I am aware, it remains a complete mystery why, for large  there
are coherent states and an apparently monostable (up to a translation) attractor,
whereas for low , there are many multistable, incoherent attractors. The argument
from Section 5.3 describes the development of a spatial discontinuity with decreasing
 but it does not explain the multistability of incoherent states. When we wrote this
paper, we did not know about the application of group theory to cluster synchrony.
The incoherent patterns are probably subgroup clusters, and it may be possible to
say something about their stability now.
106
Chapter 6: Unanswered questions
The majority of this thesis was devoted to describing what we did and why.
This Chapter focuses on what we didn’t do. The goal of this chapter is to briefly
mention some questions that became apparent during our investigation that might
lead to interesting future work.
6.1 Entropy and the photon-counting feedback loop
Time-delayed feedback loops have an infinite-dimensional phase space, and
support high-dimensional dynamics. Generally speaking, the larger the time delay
is in comparison to the filter time constants, the higher the dimensionality of the
dynamics. To obtain the results presented in Chapters 2 and 3, I deliberately
chose long filter time constants in order to observe low-dimensional dynamics. Time
series analysis techniques that rely on phase-space reconstruction, like the Cohen-
Procaccia procedure, tend not to perform well on high dimensional data, because
one needs an exceptionally large amount of samples to fill phase space sufficiently
to resolve the shape of the attractor.
From the point of view of random number generation, high-dimensional dy-
namics may offer additional security compared to low dimensional systems. An
107
hypothetical attacker who wants to predict the future state of a high-dimensional
system might try to perform a phase-space embedding. Or, in other words, they
might try to count the number of occurrences of pattern of length d, to determine
which patterns are most common. The high-dimensionality of the system will pose
an obstacle to these methods.
On the other hand, one might be able to predict the dynamics of one chaos-
based random number generator by synchronizing a second chaotic system to it.
This leads to several interesting questions.
The first question is, how does the synchronization error between two photon-
counting feedback loops depend on the photon rate λ0? There may be an information-
theoretic argument that gives some insight into the synchronization error as a
function of λ0. The synchronization could be measured using several information-
theoretic quantities including mutual information, the Kullback-Liebler divergence,
or the Jensen-Shannon divergence. It’s not clear which would be best.
The synchronization of two chaotic feedback loops could be used to construct
an upper bound for the entropy rate of a random number generator. Consider two
unidrectionally coupled feedback loops. The output of the first, which is our random
number generator, will be called Xn. The second, which will try to synchronize to
Xn, and therefore predict it, will be called Yn. Yn will be the output of a feedback
loop driven by all of the previous random numbers Xn−1, Xn−2, .... In other words,
the second loop, Y is a model that one can use to predict the next sample, Xn. If
the two loops are partially synchronized, then the conditional entropy H(Xn|Yn) <
H(Xn). In general, H(Xn|Yn) ≥ H(Xn|Xn−1, Xn−2, . . . ). So it is clearly true that
108
h ≤ (1/τ)H(X|Y ). However, because phase space embeddings of high-dimensional
signals don’t work very well, the conditional entropy might give a tighter bound on
the entropy rate than the Cohen-Procaccia method.
Another issue that could be addressed by synchronization is the effect of noise
on the unpredictability of the dynamics. When we started the work on the photon-
counting feedback loop, one of our ideas was provide an explicit demonstration that
the sensitivity to initial conditions present in chaos amplifies microscopic noise. We
still have not done this. However, it is possible to measure finite-time Lyapunov
exponents using synchronization between two loops [141]. It should also be possible
to perform an analogous experiment with the photon-counting system. This would
be an explicit measurement of how the uncertainty in the state of this system grows
in time. One might expect diffusive growth (∼
√
t) of the distance in phase space
between two trajectories in the shot-noise dominated regime, and exponential growth
in the high-photon rate limit.
6.2 Synchronization in networks of coupled systems
In the study of cluster synchrony, we have so far considered both Laplacian
coupling and adjacency matrix coupling. We have only done a little bit of work on
networks with directed links between nodes. It remains an open question whether
any interesting new phenomena are present in the case of directed coupling that are
not present with bidirectional coupling.
In my opinion, the biggest unresolved question regarding both chimeras and
109
cluster synchrony is whether the two phenomena are related. Chimera states are,
after all, situations in which a large synchronized cluster forms, while many smaller
clusters remain unsynchronized.
So far, we have largely avoided the issue of quantifying the number of possible
cluster-synchronous states that a network can have. We did count the number of
clusters and subgroups for a few different ensembles of random networks [142], but
so far we have not tried to quantify the number of subgroup clusters and Laplacian
merged clusters. This quantification would depend on the type of network. Scale-
free networks would have different statistics from Erdo˝s-Re´nyi, for example. All we
know is that there are probably a lot of different possible cluster synchronization
patterns possible in a given network. Some kind of quantification of these patterns
might go a long way towards understanding chimera states.
In the context of coupled-map lattice chimeras, the multistability of chimera
states and incoherent states needs to be quantified. Ref [139] argued that the mul-
tistability of the incoherent region scales exponentially with the system size. As far
as I am aware, it remains a complete mystery why, for large  there are coherent
states with an apparently monostable (up to a translation) attractor, whereas for
low , there are exponentially many multistable, incoherent attractors. The argu-
ment from Section 5.3 describes the development of a spatial discontinuity with
decreasing  but it does not explain the multistability of incoherent states. When
we wrote [32], we did not know about the application of group theory to cluster syn-
chrony. The incoherent patterns may be subgroup clusters, and it may be possible
to say something about their multiplicity and stability now.
110
Another question regarding coupled map lattice chimeras is why coherent pro-
files with K=2,3,4... are seen in the one-dimensional system, but not in the two
dimensional system. This difference can be seen in comparing the parameter-space
maps in Figure 5.2(a) and Figure 5.2(a). It is easy to see from inspection of Equa-
tions (5.2) and (5.3) that any solution of Equations (5.2) is automatically also a
solution of Equation (5.3), if one assumes a spatial profile that is homogeneous in
one direction. So, there are three possibilities: either wavenumbers higher than 1
are not stable in two dimensions, the basin of attraction of the higher wave numbers
is so small that they are never seen, or the transient times to these states are so
long that they are never seen.
On a more global scale, the the biggest question is how (and if) cluster syn-
chronization and chimera states manifest in natural and technological systems. We
now know that chimeras are possible in experiments. It remains an open question
whether chimeras exist in natural contexts. Another interesting potential connection
is that the symmetry-based analysis presented here might be relevant to predicting
instabilities in the synchronization of the power grid.
111
Bibliography
[1] Henri Poincare´. Lec¸ons de me´canique ce´leste: professe`es a´ la Sorbonne.
The´orie ge´ne´rale des perturbations planetaire. T. 1, volume 1. Gauthier-
Villars, 1905.
[2] Balthasar van der Pol. Forced oscillations in a circuit with non-linear resis-
tance.(reception with reactive triode). The London, Edinburgh, and Dublin
Philosophical Magazine and Journal of Science, 3(13):65–80, 1927.
[3] M. L. Cartwright and J. E. Littlewood. On nonlinear differential equations of
second order, I. The equation y¨ + κ(1 − y2)y˙ + y = bλκ cos(λt + a), κ large.
Journal of the London Mathematical Society, 20:180–189, 1945.
[4] Edward N. Lorenz. Deterministic nonperiodic flow. Journal of the Atmospheric
Sciences, 20(2):130–141, 1963.
112
[5] TienYien Li and James A. Yorke. Period three implies chaos. American
Mathematical Monthly, 82(10):985–992, 1975.
[6] Brian R. Hunt and Edward Ott. Defining chaos. Chaos, 25(9):–, 2015.
[7] Stanislaw M. Ulam and John von Neumann. On combination of stochastic
and deterministic processes-preliminary report. In Bulletin of the American
Mathematical Society, volume 53, pages 1120–1120, 1947.
[8] M. Ranferi Gutie´rrez, M. A. Reyes, and H. C. Rosu. A note on verhulst’s
logistic equation and related logistic maps. Journal of Physics A: Mathematical
and Theoretical, 43(20):205204, 2010.
[9] Per Martin-Lo¨f. The definition of random sequences. Information and Control,
9(6):602–619, 1966.
[10] Guido Boffetta, Massimo Cencini, Massimo Falcioni, and Angelo Vulpiani.
Predictability: a way to characterize complexity. Physics Reports, 356(6):367–
474, 2002.
[11] Atsushi Uchida, Kazuya Amano, Masaki Inoue, Kunihito Hirano, Sunao Naito,
Hiroyuki Someya, Isao Oowada, Takayuki Kurashige, Masaru Shiki, Shigeru
Yoshimori, Yoshimura Kazuyuk, and Peter Davis. Fast physical random bit
generation with chaotic semiconductor lasers. Nature Photonics, 2(12):728–
732, 2008.
113
[12] Ryohsuke Sakuraba, Kento Iwakawa, Kazutaka Kanno, and Atsushi Uchida.
Tb/s physical random bit generation with bandwidth-enhanced chaos in three-
cascaded semiconductor lasers. Optics Express, 23(2):1470–1490, 2015.
[13] James F. Dynes, Zhiliang L. Yuan, Andrew W. Sharpe, and Andrew J. Shields.
A high speed, postprocessing free, quantum random number generator. Ap-
plied Physics Letters, 93(3):031109, 2008.
[14] Thomas Jennewein, Ulrich Achleitner, Gregor Weihs, Harald Weinfurter, and
Anton Zeilinger. A fast and compact quantum random number generator.
Review of Scientific Instruments, 71(4):1675–1680, 2000.
[15] ID Quantique. White paper: Random number generation using quantum
physics. http://www.idquantique.com/wordpress/wp-content/uploads/
quantis-whitepaper.pdf, 2010. accessed July 30, 2015.
[16] Michael Wahl, Matthias Leifgen, Michael Berlin, Tino Ro¨hlicke, Hans-Ju¨rgen
Rahn, and Oliver Benson. An ultrafast quantum random number generator
with provably bounded output bias based on photon arrival time measure-
ments. Applied Physics Letters, 98(17):171105, 2011.
[17] Laurent Larger. Complexity in electro-optic delay dynamics: modelling, de-
sign and applications. Philosophical Transactions of the Royal Society A:
Mathematical,Physical and Engineering Sciences, 371(1999), 2013.
[18] Thomas E. Murphy, Adam B. Cohen, Bhargava Ravoori, Karl R.B. Schmitt,
Anurag V. Setty, Francesco Sorrentino, Caitlin R.S. Williams, Edward Ott,
114
and Rajarshi Roy. Complex dynamics and synchronization of delayed-feedback
nonlinear oscillators. Philosophical Transactions of the Royal Society A: Math-
ematical, Physical and Engineering Sciences, 368(1911):343–366, 2010.
[19] Y. Chembo Kouomou, Pere Colet, Laurent Larger, and Nicolas Gastaud.
Chaotic Breathers in Delayed Electro-Optical Systems. Physical Review Let-
ters, 95:203903, 2005.
[20] Michael Peil, Maxime Jacquot, Yanne Kouomou Chembo, Laurent Larger,
and Thomas Erneux. Routes to chaos and multiple time scale dynamics in
broadband bandpass nonlinear delay electro-optic oscillators. Physical Review
E, 79:026208, 2009.
[21] Andrew Rukhin, Juan Soto, James Nechvatal, Miles Smid, Elaine Barker, Ste-
fan Leigh, Mark Levenson, Mark Vangel, David Banks, Alan Heckert, et al.
NIST special publication 800-22. A statistical test suite for random and pseu-
dorandom number generators for cryptographic applications, 2001.
[22] Georges Marsaglia. DIEHARD Test suite. http://www.stat.fsu.edu/pub/
diehard/, 1998. accessed July 30, 2015.
[23] O. D’Huys, R. Vicente, T. Erneux, J. Danckaert, and I. Fischer. Synchroniza-
tion properties of network motifs: Influence of coupling delay and symmetry.
Chaos, 18:037116, 2008.
115
[24] V. Nicosia, M. Valencia, M. Chavez, A. Dı´az-Guilera, and V. Latora. Remote
synchronization reveals network symmetries and functional modules. Physical
Review Letters, 110:174102, 2013.
[25] G. Russo and J-J. E. Slotine. Symmetries, stability, and control in nonlinear
systems and networks. Physical Review E, 84:041929, 2011.
[26] C. Zhou and J. Kurths. Hierarchical synchronization in complex networks
with heterogeneous degrees. Chaos, 16(1):015104, 2006.
[27] L. M. Pecora and T. L. Carroll. Master stability functions for synchronized
coupled systems. Physical Review Letters, 80(10):2109–2112, 1998.
[28] Kuramoto, Y. . Chemical Oscillations, Waves and Turbulence. Springer-
Verlag, Berlin, 1984.
[29] Y. Kuramoto and D. Battogtokh. Coexistence of Coherence and Incoherence
in Nonlocally Coupled Phase Oscillators. Nonlinear Phenomena in Complex
Systems, 5(4):380–385, 2002.
[30] Daniel M. Abrams and S. H. Strogatz. Chimera states for coupled oscillators.
Physical Review Letters, 93(17):174102, 2004.
[31] Shin-ichiro Shima and Yoshiki Kuramoto. Rotating spiral waves with
phase-randomized core in nonlocally coupled oscillators. Physical Review E,
69(3):036213, 2004.
116
[32] Aaron M. Hagerstrom, Thomas E. Murphy, Rajarshi Roy, Philipp Ho¨vel, Iryna
Omelchenko, and Eckehard Scho¨ll. Experimental observation of chimeras in
coupled-map lattices. Nature Physics, 8(9):658–661, 2012.
[33] Mark R. Tinsley, Simbarashe Nkomo, and Kenneth Showalter. Chimera and
phase-cluster states in populations of coupled chemical oscillators. Nature
Physics, 8(9):662–665, 2012.
[34] Lennart Schmidt and Katharina Krischer. Clustering as a prerequisite
for chimera states in globally coupled systems. Physical Review Letters,
114:034101, 2015.
[35] Lennart Schmidt, Konrad Scho¨nleber, Katharina Krischer, and Vladimir
Garc´ıa-Morales. Coexistence of synchrony and incoherence in oscillatory me-
dia under nonlinear global coupling. Chaos, 24(1):013102, 2014.
[36] Erik Andreas Martens, Shashi Thutupalli, Antoine Fourrie`re, and Oskar Hal-
latschek. Chimera states in mechanical oscillator networks. Proceedings of the
National Academy of Sciences, 110(26):10563–10567, 2013.
[37] Laurent Larger, Bogdan Penkovsky, and Yuri Maistrenko. Virtual chimera
states for delayed-feedback systems. Physical Review Letters, 111(5):054103,
2013.
[38] Laurent Larger, Bogdan Penkovsky, and Yuri Maistrenko. Laser chimeras as
a paradigm for multistable patterns in complex systems. Nature Communica-
tions, 6, 2015.
117
[39] Robert McCredie May. Stability and complexity in model ecosystems, volume 6.
Princeton University Press, 2001.
[40] Leonard Mandel and Emil Wolf. Optical coherence and quantum optics. 1995.
Equation (4) is found in section 14.9.2.
[41] Y. Kouomou Chembo, Laurent Larger, Ryad Bendoula, and Pere Colet. Ef-
fects of gain and bandwidth on the multimode behavior of optoelectronic mi-
crowave oscillators. Optics Express, 16(12):9067–9072, 2008.
[42] Lennert Appeltant, Miguel Cornelles Soriano, Guy Van der Sande, Jan Danck-
aert, Serge Massar, Joni Dambre, Benjamin Schrauwen, Claudio R Mirasso,
and Ingo Fischer. Information processing using a single dynamical node as
complex system. Nature Communications, 2:468, 2011.
[43] Apostolos Argyris, Dimitris Syvridis, Laurent Larger, Valerio Annovazzi-Lodi,
Pere Colet, Ingo Fischer, Jordi Garc´ıa-Ojalvo, Claudio R Mirasso, Luis Pes-
quera, and K Alan Shore. Chaos-based communications at high bit rates using
commercial fibre-optic links. Nature, 438(7066):343–346, 2005.
[44] Bhargava Ravoori, Adam B. Cohen, Jie Sun, Adilson E. Motter, Thomas E.
Murphy, and Rajarshi Roy. Robustness of optimal synchronization in real
networks. Physical Review Letters, 107(3):034102, 2011.
[45] Caitlin R. S. Williams, Francesco Sorrentino, Thomas E. Murphy, and Ra-
jarshi Roy. Synchronization states and multistability in a ring of periodic
118
oscillators: Experimentally variable coupling delays. Chaos, 23(4):043117,
2013.
[46] Caitlin R. S. Williams, Thomas E. Murphy, Rajarshi Roy, Francesco Sor-
rentino, Thomas Dahms, and Eckehard Scho¨ll. Experimental Observations of
Group Synchrony in a System of Chaotic Optoelectronic Oscillators. Physical
Review Letters, 110:064104, 2013.
[47] Adam B. Cohen. Synchronization and prediction of chaotic dynamics on net-
works of optoelectronic oscillators, 2011.
[48] Bhargava Ravoori. Synchronization of chaotic optoelectronic oscillators:
Adaptive techniques and the design of optimal networks, 2011.
[49] Caitlin Rose Sanford Williams. Optoelectronic experiments on random bit
generators and coupled dynamical systems, 2013.
[50] J. Doyne Farmer. Chaotic attractors of an infinite-dimensional dynamical
system. Physica D: Nonlinear Phenomena, 4(3):366–393, 1982.
[51] Aaron Morgan Hagerstrom, Thomas Edward Murphy, and Rajarshi Roy. Har-
vesting entropy and quantifying the transition from noise to chaos in a photon-
counting feedback loop. Proceedings of the National Academy of Sciences,
112(30):9258–9263, 2015.
[52] Elaine Barker and John Kelsey. NIST DRAFT Special Publication 800-90b
recommendation for the entropy sources used for random bit generation. 2012.
119
[53] Nicholas Metropolis and Stanislaw Ulam. The Monte Carlo Method. Journal
of the American Statistical Association, 44(247):335–341, 1949.
[54] Arjen Lenstra, James P Hughes, Maxime Augier, Joppe Willem Bos, Thorsten
Kleinjung, and Christophe Wachter. Ron was wrong, Whit is right. Technical
report, IACR, 2012.
[55] Zvi Gutterman, Benny Pinkas, and Tzachy Reinman. Analysis of the linux
random number generator. In Security and Privacy, 2006 IEEE Symposium
on, pages 15–pp. IEEE, 2006.
[56] Mike Hamburg, Paul Kocher, and Mark E. Marson. Analysis of Intels Ivy
Bridge digital random number generator. http://www.cryptography.com/
public/pdf/Intel_TRNG_Report_20120312.pdf, 2012. accessed July 30,
2015.
[57] Taiki Yamazaki and Atsushi Uchida. Performance of random number gen-
erators using noise-based superluminescent diode and chaos-based semicon-
ductor lasers. IEEE Journal of Selected Topics in Quantum Electronics,
19(4):0600309–0600309, 2013.
[58] David P. Rosin, Damien Rontani, and Daniel J. Gauthier. Ultrafast physi-
cal generation of random numbers using hybrid Boolean networks. Physical
Review E, 87:040902, 2013.
120
[59] Reidler, I. and Aviad, Y. and Rosenbluh, M. and Kanter, I. Ultrahigh-Speed
Random Number Generation Based on a Chaotic Semiconductor Laser. Phys-
ical Review Letters, 103(2):024102, 2009.
[60] Ido Kanter, Yaara Aviad, Igor Reidler, Elad Cohen, and Michael Rosenbluh.
An optical ultrafast random bit generator. Nature Photonics, 4(1):58–61, 2009.
[61] Neus Oliver, Miguel C. Soriano, David W. Sukow, and Ingo Fischer. Dynamics
of a semiconductor laser with polarization-rotated feedback and its utilization
for random bit generation. Optics Letters, 36(23):4632–4634, 2011.
[62] Neus Oliver, Miguel Cornelles Soriano, David W. Sukow, and Ingo Fischer.
Fast random bit generation using a chaotic laser: approaching the information
theoretic limit. IEEE Journal of Quantum Electronics, 49(11):910–918, 2013.
[63] Martin Virte, Emeric Mercier, Hugo Thienpont, Krassimir Panajotov, and
Marc Sciamanna. Physical random bit generation from chaotic solitary laser
diode. Optics Express, 22(14):17271–17280, 2014.
[64] Satoshi Sunada, Takahisa Harayama, Peter Davis, Ken Tsuzuki, Ken-ichi Arai,
Kazuyuki Yoshimura, and Atsushi Uchida. Noise amplification by chaotic dy-
namics in a delayed feedback laser system and its application to nondetermin-
istic random bit generation. Chaos: An Interdisciplinary Journal of Nonlinear
Science, 22(4):047513, 2012.
[65] Takuya Mikami, Kazutaka Kanno, Kota Aoyama, Atsushi Uchida, Tohru
Ikeguchi, Takahisa Harayama, Satoshi Sunada, Ken-ichi Arai, Kazuyuki
121
Yoshimura, and Peter Davis. Estimation of entropy rate in a fast physical
random-bit generator using a chaotic semiconductor laser with intrinsic noise.
Physical Review E, 85(1):016211, 2012.
[66] Takahisa Harayama, Satoshi Sunada, Kazuyuki Yoshimura, Jun Muramatsu,
Ken-ichi Arai, Atsushi Uchida, and Peter Davis. Theory of fast nondetermin-
istic physical random-bit generation with chaotic lasers. Physical Review E,
85:046215, 2012.
[67] Xiaowen Li, Adam B. Cohen, Thomas E. Murphy, and Rajarshi Roy. Scalable
parallel physical random number generator based on a superluminescent LED.
Optics Letters, 36(6):1020–1022, 2011.
[68] Caitlin R. S. Williams, Julia C. Salevan, Xiaowen Li, Rajarshi Roy, and
Thomas E. Murphy. Fast physical random number generator using ampli-
fied spontaneous emission. Optics Express, 18(23):23584–23597, 2010.
[69] Mario Stipcˇevic´ and Daniel J. Gauthier. Precise Monte Carlo simulation
of single-photon detectors. In SPIE Defense, Security, and Sensing, pages
87270K–87270K. International Society for Optics and Photonics, 2013.
[70] Christian Gabriel, Christoffer Wittmann, Denis Sych, Ruifang Dong, Wolf-
gang Mauerer, Ulrik L. Andersen, Christoph Marquardt, and Gerd Leuchs.
A generator for unique quantum random numbers based on vacuum states.
Nature Photonics, 4(10):711–715, 2010.
122
[71] Antoine Be´rut, Artak Arakelyan, Artyom Petrosyan, Sergio Ciliberto, Raoul
Dillenschneider, and Eric Lutz. Experimental verification of Landauer’s prin-
ciple linking information and thermodynamics. Nature, 483(7388):187–189,
2012.
[72] James Sethna. Statistical mechanics: entropy, order parameters, and complex-
ity, volume 14. Oxford University Press, 2006.
[73] Hugo Touchette. The large deviation approach to statistical mechanics.
Physics Reports, 478(1):1–69, 2009.
[74] Christopher Jarzynski. Diverse phenomena, common themes. Nature Physics,
11(2):105–107, 2015.
[75] E´dgar Rolda´n and Juan M. R. Parrondo. Estimating Dissipation from Single
Stationary Trajectories. Physical Review Letters, 105:150607, 2010.
[76] M. Cencini, M. Falcioni, E. Olbrich, H. Kantz, and A. Vulpiani. Chaos or
noise: Difficulties of a distinction. Physical Review E, 62(1):427, 2000.
[77] J. G. Caputo and P. Atten. Metric entropy: An experimental means for
characterizing and quantifying chaos. Physical Review A, 35:1311–1316, 1987.
[78] C. P. Dettmann, E. G. D. Cohen, and Van Beijeren H. Statistical mechanics:
Microscopic chaos from Brownian motion? Nature, 401(6756):875–875, 1999.
123
[79] P. Gaspard, M. E. Briggs, M. K. Francis, J. V. Sengers, R. W. Gammon, J. R.
Dorfman, and R. V. Calabrese. Experimental evidence for microscopic chaos.
Nature, 394(6696):865–868, 1998.
[80] Pierre Gaspard and Xiao-Jing Wang. Noise, chaos, and (ε, τ)-entropy per unit
time. Physics Reports, 235(6):291–343, 1993.
[81] David Andrieux, Pierre Gaspard, Sergio Ciliberto, Nicolas Garnier, Sylvain
Joubaud, and Artyom Petrosyan. Entropy production and time asymmetry
in nonequilibrium fluctuations. Physical Review Letters, 98(15):150601, 2007.
[82] D. Andrieux, P. Gaspard, S. Ciliberto, N. Garnier, S. Joubaud, and A. Pet-
rosyan. Thermodynamic time asymmetry in non-equilibrium fluctuations.
Journal of Statistical Mechanics: Theory and Experiment, 2008(01):P01002,
2008.
[83] Thomas M. Cover and Joy A. Thomas. Elements of information theory. John
Wiley & Sons, 2012.
[84] Dean Prichard and James Theiler. Generalized redundancies for time series
analysis. Physica D: Nonlinear Phenomena, 84(3):476–493, 1995.
[85] Aviad Cohen and Itamar Procaccia. Computing the Kolmogorov entropy
from time signals of dissipative and conservative dynamical systems. Physical
Review A, 31:1872–1882, 1985.
[86] James Theiler. Estimating fractal dimension. Journal of the Optical Society
of America A, 7(6):1055–1073, 1990.
124
[87] Yakov Borisovich Pesin. Characteristic Lyapunov exponents and smooth er-
godic theory. Russian Mathematical Surveys, 32(4):55–114, 1977.
[88] J.-P. Eckmann and D. Ruelle. Ergodic theory of chaos and strange attractors.
Reviews of Modern Physics, 57(3):617, 1985.
[89] Giancarlo Benettin, Luigi Galgani, and Jean-Marie Strelcyn. Kolmogorov
entropy and numerical experiments. Physical Review A, 14(6):2338, 1976.
[90] Karlheinz Geist, Ulrich Parlitz, and Werner Lauterborn. Comparison of dif-
ferent methods for computing Lyapunov exponents. Progress of Theoretical
Physics, 83(5):875–893, 1990.
[91] Peter Grassberger and Itamar Procaccia. Measuring the strangeness of strange
attractors. Physica D, 9:189–208, 1983.
[92] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa. Spontaneous syn-
chrony in power-grid networks. Nature Physics, 9(3):191–197, 2013.
[93] A.-L. Do, Johannes Ho¨fener, and Thilo Gross. Engineering mesoscale
structures with distinct dynamical implications. New Journal of Physics,
14(11):115022, 2012.
[94] T. Dahms, J. Lehnert, and E. Scho¨ll. Cluster and group synchronization in
delay-coupled networks. Physical Review E, 86:016202, 2012.
125
[95] Chenbo Fu, Zhigang Deng, Liang Huang, and Xingang Wang. Topological
control of synchronous patterns in systems of networked chaotic oscillators.
Physical Review E, 87(3):032909, 2013.
[96] I. Kanter, M. Zigzag, A. Englert, F. Geissler, and W. Kinzel. Synchronization
of unidirectional time delay chaotic networks and the greatest common divisor.
Europhysics Letters, 93:6003, 2011.
[97] D. P. Rosin, D. Rontani, D. J. Gauthier, and E. Scho¨ll. Control of synchro-
nization patterns in neural-like boolean networks. Physical Review Letters,
110(10):104102, 2013.
[98] F. Sorrentino and E. Ott. Network synchronization of groups. Physical Review
E, 76:056114, 2007.
[99] Vladimir N. Belykh, Grigory V. Osipov, Valentin S. Petrov, Johan A. K.
Suykens, and Joos Vandewalle. Cluster synchronization in oscillatory net-
works. Chaos, 18:037106, 2008.
[100] M. Golubitsky, I. Stewart, and D.G. Schaeffer. Singularities and groups in
bifurcation theory, volume II. Springer-Verlag, New York, NY, 1985.
[101] M. Golubitsky and I. Stewart. The Symmetry Perspective: From Equilibrium
to Chaos in Phase Space and Physical Space. Berkha¨user-Verlag, Basel, 2002.
[102] Kevin Judd. Networked dynamical systems with linear coupling: Synchroni-
sation patterns, coherence and other behaviours. Chaos, 23:043112, 2013.
126
[103] M. Newman. Networks, An Introduction. Oxford University Press, 2011. ,
Chap.18.
[104] Kenneth S. Fink, Gregg Johnson, Tom Carroll, Doug Mar, and Lou Pec-
ora. Three-oscillator systems as universal probes of coupled oscillator stability.
Physical Review E, 61(5):5080–5090, 2000.
[105] D. Irving and F. Sorrentino. Synchronization of dynamical hypernetworks:
dimensionality reduction through simultaneous block-diagonalization of ma-
trices. Physical Review E, 86:056102, 2012.
[106] Alex Arenas, Jurgen Kurths Albert Dı´az-Guilera, Yamir Moreno, and Chang-
song Zhou. Synchronization in complex networks. Physics Reports, 469:93–
153, 2008.
[107] M. Tinkham. Group Theory and Quantum Mechanics. McGraw-Hill, New
York, NY, 1964.
[108] William Stein. SAGE: Software for Algebra and Geometry Experimentation.
http://www.sagemath.org/sage/, http://sage.scipy.org/, 2013.
[109] The GAP Group. GAP: Groups, Algorithms, and Programming, Version 4.4.
http://www.gap-system.org, 2005.
[110] F. Vallentin. Symmetry in semidefinite programs. Linear Algebra Appl.,
430(1):360–369, 2009.
127
[111] R. Goodman and N.R. Wallach. Representations and Invariants of the Clas-
sical Groups. Cambridge University Press, 1998.
[112] B. E. Sagan. The Symmetric Group. Wadsworth Brooks, Pacific Grove, CA,
1991.
[113] B. D. MacArthur and R. J. Sa´nchez-Garc´ıa. Spectral characteristics of network
redundancy. Physical Review E, 80:026117, 2009.
[114] B. D. MacArthur, R.J. Sa´nchez-Garc´ıa, and J. W. Anderson. On automor-
phism groups of networks. Discrete Applied Mathematics, 156:3525–3531,
2008.
[115] B. Dionne, M. Golubitsky, and I. Stewart. Coupled cells with internal sym-
metry: I. wreath products. Nonlinearity, 9:559–574, 1996.
[116] Kunihiko Kaneko. Clustering, coding, switching, hierarchical ordering, and
control in a network of chaotic elements. Physica D: Nonlinear Phenomena,
41(2):137–172, 1990.
[117] Francesco Sorrentino, Louis M. Pecora, Aaron M. Hagerstrom, Thomas E.
Murphy, and Rajarshi Roy. Cluster Synchronization in Laplacian Coupled
Systems: Beyond Symmetries. Submitted to Physical Review Letters (arXiv
preprint arXiv:1507.04381).
[118] Edward Ott and Thomas M. Antonsen. Low dimensional behavior of large
systems of globally coupled oscillators. Chaos, 18(3):037113, 2008.
128
[119] Daniel M. Abrams, Rennie Mirollo, S. H. Strogatz, and Daniel A. Wiley. Solv-
able model for chimera states of coupled oscillators. Physical Review Letters,
101(8):084103, 2008.
[120] Erik A. Martens, Carlo. R. Laing, and Steven. H. Strogatz. Solvable model of
spiral wave chimeras. Physical Review Letters, 104(4):044101, 2010.
[121] Adilson E. Motter. Nonlinear dynamics: Spontaneous synchrony breaking.
Nature Physics, 6(3):164–165, 2010.
[122] C. R. Laing. The dynamics of chimera states in heterogeneous Kuramoto
networks. Physica D, 238(16):1569–1588, 2009.
[123] C. R. Laing. Chimeras in networks of planar oscillators. Physical Review E,
81(6):066221, 2010.
[124] C. R. Laing. Fronts and bumps in spatially extended Kuramoto networks.
Physica D, 240(24):1960–1971, 2011.
[125] Erik A. Martens. Bistable chimera attractors on a triangular network of os-
cillator populations. Physical Review E, 82:016216, 2010.
[126] Erik A. Martens. Chimeras in a network of three oscillator populations with
varying network topology. Chaos: An Interdisciplinary Journal of Nonlinear
Science, 20(4):043122, 2010.
[127] Matthias Wolfrum and Oleh E. Omel’chenko. Chimera states are chaotic
transients. Physical Review E, 84:015201, 2011.
129
[128] Gautam C. Sethia, Abhijit Sen, and Fatihcan M. Atay. Clustered chimera
states in delay-coupled oscillator systems. Physical Review Letters, 100:144102,
2008.
[129] Abhijit Sen, Ramana Dodla, George L. Johnston, and Gautam C. Sethia.
Understanding complex systems. (ed.) Atay, F. M. pages pp. 1–41. Springer-
Verlag, Berlin, 2010.
[130] Irene Waller and Raymond Kapral. Spatial and temporal structure in systems
of coupled nonlinear oscillators. Physical Review A, 30:2047–2055, 1984.
[131] K. Kaneko. Overview of coupled map lattices. Chaos, 2:279 –282, 1992.
[132] K. Kaneko. Theory and applications of coupled map lattices. John Wiley &
Sons, Chichester, 1993.
[133] K. Kaneko. Period-doubling of kink-antikink patterns, quasiperiodicity in
anti-ferro-like structures and spatial intermittency in coupled logistic lattice.
Progress of Theoretical Physics, 72(3):480–486, 1984.
[134] Arkady S. Pikovsky, Michael Rosenblum, and Ju¨rgen Kurths. Synchronization:
A universal concept in nonlinear sciences, volume 12. Cambridge Univerity
Press, Cambridge, 2003.
[135] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, and C. S. Zhou. The
synchronization of chaotic systems. Physics Reports, 366:1–101, 2002.
130
[136] Carlo R. Laing. Chimeras in networks with purely local coupling. arXiv
preprint arXiv:1506.05871, 2015.
[137] Gautam C. Sethia and Abhijit Sen. Chimera states: the existence criteria
revisited. Physical Review Letters, 112(14):144101, 2014.
[138] Mark J. Panaggio and Daniel M. Abrams. Chimera states on a flat torus.
Physical Review Letters, 110(9):094102, 2013.
[139] Iryna Omelchenko, Yuri L. Maistrenko, Philipp Ho¨vel, and Eckehard Scho¨ll.
Loss of coherence in dynamical networks: spatial chaos and chimera states.
Physical Review Letters, 106:234102, 2011.
[140] Iryna Omelchenko, Bruno Riemenschneider, Philipp Ho¨vel, Yuri L.
Maistrenko, and Eckehard Scho¨ll. Transition from spatial coherence to in-
coherence in coupled chaotic systems. Physical Review E, 85:026212, 2012.
[141] Adam B. Cohen, Bhargava Ravoori, Thomas E. Murphy, and Rajarshi Roy.
Using synchronization for prediction of high-dimensional chaotic dynamics.
Physical Review Letters, 101(15):154102, 2008.
[142] Louis M. Pecora, Francesco Sorrentino, Aaron M. Hagerstrom, Thomas E.
Murphy, and Rajarshi Roy. Cluster synchronization and isolated desynchro-
nization in complex networks with symmetries. Nature Communications, 5,
2014.
131