Stochastic and Collective Properties of

Stochastic and Collective Properties of Nonlinear Oscillators - CaltechTHESIS
CaltechTHESIS
A Caltech Library Service
About
Browse
Deposit an Item
Instructions for Students
Stochastic and Collective Properties of Nonlinear Oscillators
Citation
Kogan, Oleg Boris
(2009)
Stochastic and Collective Properties of Nonlinear Oscillators.
Dissertation (Ph.D.), California Institute of Technology.
doi:10.7907/93R8-TJ70.
Abstract
Two systems of nonlinear oscillators are considered: (a) a single periodically driven nonlinear oscillator interacting with a heat bath, which may operate in the regime of bistability or monostability, and (b) a one-dimensional chain of self-sustaining phase oscillators with nearest-neighbor interaction.
For a single oscillator we analyze the scaling crossovers in the thermal activation barrier between the two stable states. The rate of metastable decay in nonequilibrium systems is expected to display scaling behavior: the logarithm of the decay rate should scale as a power of the distance to a bifurcation point where the metastable state disappears. We establish the range where different scaling behavior is displayed and show how the crossover between different types of scaling occurs. Using the instanton method, we map numerically the entire parameter range of bistability and find the regions where the scaling exponents are 1 or 3/2, depending on the damping. The exponent 3/2 is found to extend much further from the bifurcation then where it would be expected to hold as a result of an overdamped soft mode. Additionally, we uncover a new scaling behavior with exponent of ≈1.3 that extends beyond the close vicinity of the bifurcation point.
We also study the pattern of fluctuational trajectories in the monostable regime. For nonequilibrium systems, fluctuational and relaxational trajectories are not simply related by time-reversibility, as is the case in thermal equilibrium. One of the consequences of this is the onset of singularities in the pattern of fluctuational trajectories, where most probable paths to neighboring states are far away from each other. This also creates nonsmoothness in the probability distribution of the system in its phase space. We discover that the pattern of optimal paths in equilibrium systems is fragile with respect to the driving strength F, and investigate how the singularities occur as the system is driven away from equilibrium. As the strength of the driving F approaches zero, the cusp of the spiral caustic system recedes to larger radius R and the angle of the cusp also decreases. The dependence of R on F displays two scaling laws with crossovers, where the scaling exponents depend on the damping.
For the one-dimensional chain of nearest-neighbor coupled phase oscillators, we develop a renormalization group method to investigate synchronization clusters. We apply it numerically to Lorentzian distributions of intrinsic frequencies and couplings and investigate the statistics of the resultant cluster sizes and frequencies. We find that the distributions of sizes of frequency clusters are exponential, with a characteristic length. The dependence of this length upon parameters of these Lorentzian distributions develops an asymptotic power law with an exponent of 0.48 ± 0.02. The findings obtained with the renormalization group are compared with numerical simulations of the equations of motion of the chain, with an excellent agreement in all the aforementioned quantities.
Item Type:
Thesis (Dissertation (Ph.D.))
Subject Keywords:
collective; disorder; fluctuations; nonequilibrium; nonlinear; randomness; rare events
Degree Grantor:
California Institute of Technology
Division:
Engineering and Applied Science
Major Option:
Materials Science
Thesis Availability:
Public (worldwide access)
Research Advisor(s):
Cross, Michael Clifford
Thesis Committee:
Cross, Michael Clifford (chair)
Fultz, Brent T.
Refael, Gil
Corngold, Noel Robert
Johnson, William Lewis
Defense Date:
16 December 2008
Record Number:
CaltechETD:etd-06012009-145134
Persistent URL:
DOI:
10.7907/93R8-TJ70
Default Usage Policy:
No commercial reproduction, distribution, display or performance rights in this work are provided.
ID Code:
2370
Collection:
CaltechTHESIS
Deposited By:
Imported from ETD-db
Deposited On:
02 Jun 2009
Last Modified:
26 Nov 2019 19:13
Thesis Files
Preview
PDF
- Final Version
See Usage Policy.
22MB
Repository Staff Only:
item control page
CaltechTHESIS is powered by
EPrints 3.3
which is developed by the
School of Electronics and Computer Science
at the University of Southampton.
More information and software credits
Stochastic and Collective Properties of Nonlinear
Oscillators

Thesis by

Oleg Kogan
In Partial Fulfillment of the Requirements
for the Degree of
Doctor of Philosophy

California Institute of Technology
Pasadena, California

2009
(Defended December 16, 2008)

c 2009
Oleg Kogan

iii

Dedicated to the memory of my mother.

Acknowledgements
I would like to thank my advisor and teacher Michael Cross for his patience and
countless chats about physics through all these years. I always found these conversations to be very insightful and would like to hope that they had an influence on
my own development as a scientist. I always aspire to attempt to see problems in a
similarly intuitive and broad point of view. I also thank Mark Dykman for his “unofficial mentorship” and support throughout my many visits to MSU and for many
also insightful, long phone conversations that we have had about physics. I also
thank Michael Khasin for useful discussions. Part II of this work would not have
been possible without the clever renormalization ideas of Gil Refael and his general
expert-level understanding of renormalization group and critical phenomena. I feel
lucky to have come across these three very strong physicists and look forward to
potential future collaborations. The work described in part II would also have been
much more difficult without the input from Jeffrey Rogers, especially in the numerical
simulations. On the personal, but not any less important side, I thank my father for
making it possible for me to be able to go to college and then to graduate school in
the first place, while enduring and overcoming the hardships of assimilation in a new
country that can only be understood by few. I also thank Sima for her kindness and
support. I warmly acknowledge Josephine’s presence through these years. Finally, I
hold my officemates Janet Scheel and Tony Lee accountable for bringing nonscientific
discussions into our office life.

Abstract
Two systems of nonlinear oscillators are considered: (a) a single periodically driven
nonlinear oscillator interacting with a heat bath, which may operate in the regime of
bistability or monostability, and (b) a one-dimensional chain of self-sustaining phase
oscillators with nearest-neighbor interaction.
For a single oscillator we analyze the scaling crossovers in the thermal activation
barrier between the two stable states. The rate of metastable decay in nonequilibrium
systems is expected to display scaling behavior: the logarithm of the decay rate should
scale as a power of the distance to a bifurcation point where the metastable state
disappears. We establish the range where different scaling behavior is displayed and
show how the crossover between different types of scaling occurs. We map the entire
parameter range of bistability and find the regions where the scaling exponents are
1 or 3/2, depending on the damping. The exponent 3/2 is found to extend much
further from the bifurcation then where it would be expected to hold as a result of
an overdamped soft mode. Additionally, we uncover a new scaling behavior with
exponent of ≈ 1.3 that extends beyond the close vicinity of the bifurcation point.
We also study the pattern of fluctuational trajectories in the monostable regime.
For nonequilibrium systems, fluctuational and relaxational trajectories are not simply
related by time-reversibility, as is the case in thermal equilibrium. One of the consequences of this is the onset of singularities in the pattern of fluctuational trajectories,
where most probable paths to neighboring states are far away from each other. This
also creates nonsmoothness in the probability distribution of the system in its phase
space. We discover that the pattern of optimal paths in equilibrium systems is fragile
with respect to the driving strength F, and investigate how the singularities occur

vi
as the system is driven away from equilibrium. As the strength of the driving F
approaches zero, the cusp of the spiral caustic system recedes to larger radius R and
the angle of the cusp also decreases. The dependence of R on F displays two scaling
laws with crossovers, where the scaling exponents depend on the damping.
For the one-dimensional chain of nearest-neighbor coupled phase oscillators, we
develop a renormalization group method to investigate synchronization clusters. We
apply it numerically to Lorentzian distributions of intrinsic frequencies and couplings
and investigate the statistics of the resultant cluster sizes and frequencies. We find
that the distributions of sizes of frequency clusters are exponential, with a characteristic length. The dependence of this length upon parameters of these Lorentzian
distributions develops an asymptotic power law with an exponent of 0.48 ± 0.02. The
findings obtained with the renormalization group are compared with numerical simulations of the equations of motion of the chain, with an excellent agreement in all
the aforementioned quantities.

vii

Contents
Acknowledgements

Abstract

1 Overview

Aspects of Thermal Activation of Nonlinear Oscillators

out of Equilibrium

2 Driven Nonlinear Oscillator Interacting with a Heat Bath

2.1

2.2

2.3

From Langevin Equation to Lagrangian Manifolds . . . . . . . . . . .

2.1.1

Slow Envelope Dynamics . . . . . . . . . . . . . . . . . . . . .

2.1.2

Fokker-Planck equation and the Path Integral solution . . . .

2.1.3

Auxiliary Classical Mechanics . . . . . . . . . . . . . . . . . .

2.1.4

Escape Rates . . . . . . . . . . . . . . . . . . . . . . . . . . .

Analytical Limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2.1

Close to Bifurcations . . . . . . . . . . . . . . . . . . . . . . .

2.2.2

Low Effective Damping Limit . . . . . . . . . . . . . . . . . .

Power Law Scalings of Escape Barriers . . . . . . . . . . . . . . . . .

3 Scaling Crossovers

3.1

Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2.1

Up → Down Transitions . . . . . . . . . . . . . . . . . . . . .

viii

3.2.2

3.2.1.1

Numerical Results . . . . . . . . . . . . . . . . . . .

3.2.1.2

Discussion of the ξ = 3/2 to ξ = 1 Crossover . . . .

Down → Up Transitions . . . . . . . . . . . . . . . . . . . . .

4 Consequence of Breaking of Detailed Balance

4.1

Detailed Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2

Fragility of Detailed Balance and Formation of Caustics in the Monos-

53
54

table Regime of the Oscillator . . . . . . . . . . . . . . . . . . . . . .

4.2.1

Exact Solution at Zero F

. . . . . . . . . . . . . . . . . . . .

4.2.2

An Aside . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2.3

Phenomenology of Caustics in the Monostable Regime . . . .

4.2.4

Analytical Predictions . . . . . . . . . . . . . . . . . . . . . .

4.3

Physics Without the Slowing-down Region . . . . . . . . . . . . . . .

4.4

Perturbation Theory for r(p) . . . . . . . . . . . . . . . . . . . . . .

Renormalization Group Method for Predicting Frequency

Clusters in a Chain of Nearest-Neighbor Phase Oscillators 87
5 Synchronization: Background

6 Oscillator Chain RG

6.1

Idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.2

Strong-Coupling Step . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.3

Crazy Oscillator Step . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.3.1
6.4

Comments on the Crazy Oscillator Step . . . . . . . . . . . . 109

Renormalization Group Implementation . . . . . . . . . . . . . . . . 113
6.4.1

Strain Check . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

7 Results

118

7.1

Testing the validity of the numerical RG . . . . . . . . . . . . . . . . 118
7.1.1

RG Analysis of Lorentzian Chains . . . . . . . . . . . . . . . . 119

7.2

7.1.2

Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

7.1.3

Comparison of RG and Numerics . . . . . . . . . . . . . . . . 120

Physics of the Unsynchronized Phase . . . . . . . . . . . . . . . . . . 122
7.2.1

Cluster Size Distribution . . . . . . . . . . . . . . . . . . . . . 122

7.2.2

Dependence of ξ on Disorder Strength . . . . . . . . . . . . . 125

7.2.3

Approximate Analytical Argument . . . . . . . . . . . . . . . 128

7.2.4

Cluster Frequency Distributions . . . . . . . . . . . . . . . . . 129

8 Concluding Remarks

133

A Relating Beam Theory to Duffing Equation

137

B Properties of Noise in the Amplitude Equations

142

C Features in the Pattern of Classical Trajectories

145

C.1 Shadow regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
C.2 Caustics emanating from saddles . . . . . . . . . . . . . . . . . . . . 146
C.3 Rotation of eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . 147
D Full Locking Criterion

156

Bibliography

161

List of Figures
2.1

Amplitude of fixed points versus F for three different values of damping:
η = 0.01ηc , η = 0.5ηc , and η = ηc = √13 , left to right. . . . . . . . . . .

2.2

The dynamics in the P -Q space for a particular values of parameters
η = 0.3σc ≈ 0.173 and F = 0.3. The fixed points are shown as black
solid dots. The left-most one is the saddle, which lies on the separatrix
between two basins. The other curves shown are trajectory from the
saddle to each of the attractors (blue spiraling curves) and some other
typical trajectory (dashed curve). . . . . . . . . . . . . . . . . . . . . .

2.3

The region of bistability in the η-F parameter space is the closed region
between the FBh (top) and FBl (bottom) curves. These curves meet at
q ´
. The width of the hysteresis at
the cusp point at ηc = √13 , Fc = 27
. Outside of the bistability region there is always only one
zero η is 27
fixed point which is an attractor. . . . . . . . . . . . . . . . . . . . . .

2.4

Cartoon of the K-flow
near a saddle-node bifurcation. . . . . . . . . .

2.5

Activation barriers in the limit of η → 0 given by Eq. (2.69). The curve
that increases from left to right corresponds to Rup→down and the curve
that decreases from left to right corresponds to Rdown→up . . . . . . . .

3.1

Main plots: S ∗ (x)/η (dots) and fits (lines) for up → down transitions.
The fit to the ξ = 3/2 regime was made by analyzing ln (S ∗ /η) versus
ln x and a fit to the ξ = 1 regime was made by analyzing S ∗ /η versus
x. Inserts: S ∗ (x)/η (dots) and theory (lines) based on Eq. (2.63) for the
ξ = 3/2 regime and on the quasi-Hamiltonian theory, Eq. (2.69). . . . .

xi
3.2

S ∗ (x)/η on a linear scale. The thick black line represents the zero damping theory, Eq. (2.69), the same as one of the curves in Fig. 2.5. The
other curves represent numerical calculations of S ∗ : squares for η = 0.03
and circles for η = 0.06. Note the linear regime at finite η is delayed
until larger F, giving way to the 3/2 regime prior to it.

3.3

. . . . . . . .

(a),(b) Up → down escape trajectories for η = 0.001; (a) at 7% of the
hysteresis (i.e. x = 0.07), (b) at 0.001% of the hysteresis. (c),(d) Up →
down escape trajectories for η = 0.1; (c) at 80% of the hysteresis, (d) at
0.001% of the hysteresis. . . . . . . . . . . . . . . . . . . . . . . . . . .

3.4

Regions of different scaling behaviors for up → down transitions obtained by scanning F at fixed η for multiple values of η. The meaning
of the three scaling regions is described in the text. The dashed line is a
simple prediction of the crossover between the ξ = 1 and ξ = 3/2 regime
expressed by Eq. (3.3). . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5

S ∗ /η versus x and η up to η = 0.175 which is about 30% of the hysteresis.
Fig. 3.1 are two slices of this surface. . . . . . . . . . . . . . . . . . . .

3.6

R1 (x)/η ≡ (S ∗ (x) − R0 (x))/η plotted for η = 0.01. The straight line
indicates a fit within the range spanning from x ≈ 0.0013 to x ≈ 0.019;
the slope of the fit for this particular value of η is ≈ 2.51. . . . . . . .

3.7

Main plot: the exponent p obtained from the fit of R1 (x)/η, plotted versus η. Insert, the position of LC defined as that x at which
0.2R0 (x)/η = R1 (x)/η. . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.8

(a) S(l) and its derivative for particular parameters in region III: here
η = 0.15 and F = 10% of the hysteresis, corresponding to µ ≈ 5; the
soft mode picture is expected to hold only for µ ¿ 1. (b) The plot of
the MPEP along which S(l) in (a) is computed; dashed curve denotes
separatrix. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xii
3.9

Regions of different scaling behaviors for down → up transitions obtained by scanning F at fixed η for multiple values of η. Region I corresponds to a nonpower-law regime. Region II corresponds to the ξ ≈ 1.3
regime. Region III corresponds to the ξ = 3/2 regime. Fuzzy regions
correspond to parts of the (F, η) plane which proved very difficult to
map reliably due to the high divergence of trajectories. . . . . . . . . .

3.10

Main plots: S ∗ (x)/η (dots) and fits (lines) for down → up transitions.
Fits to the ξ = 3/2 and ξ ≈ 1.3 regimes was made by analyzing ln (S ∗ /η)
versus ln x. In the region where the simple shooting method was unreliable for hitting the saddle due to the high divergence of trajectories, a
bracketing method was used. The results from using this method are denoted by squares in plot (b). Inserts: S ∗ (x)/η (dots) and theory (lines)
based on the quasi-Hamiltonian theory, Eq. (2.69). Note: this theory
predicts the ξ = 3/2 scaling, as was first pointed out by [77]; in the
light of the recent work by [21] this was related to the nonlocality of the
bifurcation at zero η. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.11

S ∗ /η versus x and η up to η = 0.22 which is about 33% of the hysteresis.
Fig. 3.1 are two slices of this surface. . . . . . . . . . . . . . . . . . . .

4.1

The pattern of 10 escape trajectories at F = 0. Each trajectory undergoes approximately 3 revolutions after reversing its direction of rotation
at the slow-down radius of 1. . . . . . . . . . . . . . . . . . . . . . . .

4.2

The pattern of 10 escape trajectories at F = 0.005. . . . . . . . . . . .

4.3

The pattern of 10 escape trajectories at F = 0.008 and the pair of
caustics originating at the cusp. . . . . . . . . . . . . . . . . . . . . . .

4.4

Log-log plot representing the radius of the cusp, measured from the
origin versus the strength of the driving F for η = 0.4 circles and η = 0.1
diamonds. The solid and dashed lines represent fits to the two power-law
regimes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiii
4.5

Caustics systems at various F. Panels (a) and (b): F = 0.05, panels
(c) and (d): F = 0.2, and panels (e) and (f): F = .395, all at η = 0.4.
In panel (f) we also show a set of escape trajectories from one of the
attracting fixed points, one for a different trajectory parameter. The
fixed points and the separatrix are also shown, with the saddle laying
on the separatrix. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.6

Radius of the cusp versus η for a fixed F ≈ 0.0012. The asymptotic
exponent p̃ in R = 1 + Aη p̃ is ≈ 1.26. . . . . . . . . . . . . . . . . . . .

4.7

A portion of the caustic in the system based on g(Q, P ) = 41 (Q2 + P 2 + 1) ,
indicated by open circles and a sample of 5 trajectories. Note the absence of the turn-around at r = 1. . . . . . . . . . . . . . . . . . . . . .

5.1

Development of frequency clusters as the value of coupling K is increased, shown on a small section of a 1000 oscillator chain. Solid diamonds: time-averaged running frequencies of oscillators, defined similarly to Eq. (6.40) with t − t0 a large, but finite time interval; lines:
intrinsic frequencies that appear in Eqn. (5.4). . . . . . . . . . . . . . .

6.1

A small piece of the chain around the strongest coupling, labelled here
as Kn . All the neighboring Ki and all the local ωi are assumed to have
much smaller magnitude then Kn .

6.2

. . . . . . . . . . . . . . . . . . . .

A small piece of the chain around the oscillator with the largest frequency, labelled here as Ω. All the neighboring Ki and all the neighboring ωi are assumed to have much smaller magnitude then ωn . . . . . .

xiv
6.3

Top: the set of intrinsic frequencies and couplings. The frequency of
oscillator 4 has been varied. (a) ωi (Ω) for oscillator 3 (green circles),
oscillator 4 (red squares) and oscillator 5 (blue diamonds). The squareroot approximations, Eqs. (6.50),(6.51 are shown by the thick curves, the
simple formula such as Eq. (6.45) appears as the dotted curve and the
function ω̄ = Ω appears in light grey. (b) Close-up of the low-Ω region.
Also shown are the average frequencies of all three oscillators (thin solid
line), of oscillators 3 and 4 (dotted line) and of oscillators 4 and 5 (dashed
line). (c) Same as (a) when oscillators 3, 4, 5 are embedded into a 7oscillator system with frequencies and couplings specified in the top
plot. (d) ω̄4 (Ω) for three different functional forms of the interaction.
The curve closest to ω̄ = Ω grey line is the saw-tooth interaction, the one
furthest is the 4-th order polynomial interaction and the middle curve
is the sine interaction, i.e. same as in (a). . . . . . . . . . . . . . . . . 111

6.4

Triangle wave, sine (middle curve) and quartic polynomial interaction
(outer curve). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

7.1

Cluster structure of a small section of a 10000-oscillator chain. The
three columns correspond to µ = 1.25, µ = 3.75 and µ = 7.5 respectively. The top row displays 500 time-averaged running frequencies ω̄
from simulations of Eq. (6.2). The middle row zooms in on a smaller
section of 50 oscillators and compares ω̄ from simulations (solid squares)
and from the RG (open circles). The bottom row displays the difference
between the frequencies ω̄ of the RG and the simulations for the same
50 oscillators as in the middle row. . . . . . . . . . . . . . . . . . . . . 121

7.2

Number of clusters of a given size versus the size. Comparison between
simulations of Eq. (6.2) (solid squares), the RG without strain check
(open circles) and the RG with the strain check (diamonds). The curves
are spread out for clarity by multiplying each successive data set as µ
decreases by 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

xv
7.3

Number of clusters of a given size P (n), normalized by the system size
for µ = 7.5 at two system sizes: 106 oscillators (open symbols) and 107
oscillators (smaller filled symbols). (a) RG without strack check and (b)
RG with the strain check. . . . . . . . . . . . . . . . . . . . . . . . . . 124

7.4

(a) P(n) along with the fit between cluster size 7 and 35 as described
in the text. (b) The slopes of partial fits between cluster size m and
cluster size 35. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

7.5

Characteristic cluster length ξ given by the inverse of the fits to data in
Fig. 7.2. (a): Results of simulation (filled squares), the numerical RG
of Section 6.1 (open circles) and an enhanced numerical RG with strain
check (diamonds) described in Section 7.2.2. The solid line is a simple
prediction ξ(µ) = −1/ log(1 − ℘(µ)) where ℘(µ) comes from Eq. (7.8).
The dashed line is a linear fit to the simulation data (solid squares) for
µ > 0.625 and has a slope of 0.461; similar fits to the circles give the
slope 0.477 and to the diamonds the slope of 0.478. (b): Log-log plot of
the same data excluding the analytic prediction. . . . . . . . . . . . . . 126

7.6

Slopes of partial fits of ξ versus ln µ between ln µlower and ln 7.5. Solid
dots: simulation; open circles: RG without strain check; diamonds: RG
with the strain check. . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

7.7

Frequency distributions. The y-axis displays the number of clusters of a
given size m in a frequency bin of width ∆ω = 0.1 at frequency specified
on the x-axis.We display distributions for m = 1, 2, 3, 4 at µ = 1.25, 3.75
and 7.5. Thick curve: numerical simulation; thin red curve: numerical
RG without the use of the strain check. . . . . . . . . . . . . . . . . . . 130

7.8

Frequency distributions. The y-axis displays the number of clusters of a
given size m in a frequency bin of width ∆ω = 0.1 at frequency specified
on the x-axis.We display distributions for m = 1, 2, 3, 4 at µ = 1.25, 3.75
and 7.5. Thick curve: numerical simulation; thin green curve: numerical
RG with the use of the strain check. . . . . . . . . . . . . . . . . . . . 131

xvi
8.1

Fixed points and flow diagram of the oscillator chain. The short-range
oscillator chain has only two fixed points: the unstable fully synchronized fixed point at infinite interactions (or zero disorder) µ → ∞, and
the stable unsynchronized fixed point. The cross over flow between the
two points is associated with a correlation length, captured by the average cluster size, ξ(µ).

A.1

. . . . . . . . . . . . . . . . . . . . . . . . . . 134

Top: ω12 (ξ)/ξ. Notice the crossover from the the string-dominated regime
at lower ξ to the stiffness-dominated regime at larger ξ. Bottom, −γ/µ. 141

C.1

A set of 50 trajectories that escape the attracting fixed point are shown
in blue, the most likely escape trajectory (which we know must pass
through the saddle) is in red and the relaxational trajectory from the
saddle to the attracting fixed point appears in black. Note, that these
two are not just related by reversing the sign of the friction, as would be
in the case of detailed balance. The separatrix of the lower-amplitude
basin is shown in green. . . . . . . . . . . . . . . . . . . . . . . . . . . 149

C.2

Same features as in Fig. C.1, but for a narrower range of parameters,
[4.5, 5.985] and closer to the saddle. The cusp and the pair of caustics
are clearly seen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

C.3

A set of 50 trajectories in the range of parameters [5.6328, 5.63398]. The
the region shown is even smaller then in Fig. C.2. . . . . . . . . . . . . 151

C.4

Same features as in Fig. C.3, but showing the larger region of the (Q, P )
space. Some trajectories hit the separatrix very close to the saddle and
some, for a slightly different trajectory parameter, lead to a completely
different region of the (Q, P ) space. . . . . . . . . . . . . . . . . . . . . 152

C.5

A set of 50 trajectories in a very narrow window of trajectory parameters, [5.633489536, 5.63348954776]. The projection of the eigen-direction
of the most-repelling eigenvector is also shown in magenta. . . . . . . . 153

C.6

Two sets of trajectories, one for a wide parameter region and one for a
much narrower one. The saddle point is not shown. . . . . . . . . . . . 154

xvii
C.7

Pattern of escape trajectories for up → down transitions at η = 0.15.
Top to bottom: patterns of escape trajectories at F corresponding to
0.8%,0.9% and 0.989998% of the hysteresis respectively. The separatrix
(green) and the saddle point (dark blue dot) are also shown. . . . . . . 155

D.1

The functions u(m, µ/λ) and l(m, µ/λ) plotted versus µ/λ for three
values of m: 1 (solid line), 3 (short dashed curve) and 7 (long dashed
curve). In each case, the upper and lower curve corresponds, respectively
to functions u and l. Notice the logarithmic scale of the x-axis. . . . . 160

Chapter 1
Overview
Equilibrium statistical mechanics serves as an indisputably powerful and generic tool
that paved way for much progress in our understanding of collective phenomena in
nature. This success lies in contrast to what is known outside the realm of equilibrium
Hamiltonian systems. This includes driven Hamiltonian systems, or systems which
are not based on a Hamiltonian at all. Although a great variety of techniques for
nonequilibrium systems have been developed [1, 2], a unifying theory does not yet
exist. While a big fraction of frontier research in modern physics lies in the quantum
regime, understanding the laws that govern nonequilibrium systems, even the classical
ones, is also fundamentally significant. At present, the research in this direction lies
in the information-gathering stage [3]. A particularly important task is to find generic
aspects of such systems. A wide range of problems, from the motion of vortex arrays
in superconductors and condensates to the human brain, are examined with the hope
of identifying unifying concepts which may eventually pave the way to a more generic
understanding of nonequilibrium phenomena.
Nonlinear oscillators model a rich set of natural or man-made systems which can
be used to probe fundamental issues in nonequilibrium statistical mechanics. They are
also of experimental or practical value. For example, the advent of nanotechnology has
brought us to the exciting length and time scales where nonlinearity of nanooscillating
devices becomes important [4, 5] (see also appendix A), concomittantly with the rising
importance of thermal noise [4] (and references therein). It has been demonstrated
that nanoelectromechanical systems (NEMS), when built to operate in the classical

regime are accurately modelled as nonlinear oscillators (see our recent work [6] for
instance). Additionally, if the current trend in miniaturization continues, it is very
likely that quantum effects, such as tunneling, discretization of energy levels, quantum
uncertainties, and entanglement are going to become apparent in the foreseeable
future (see for example [7] - [11] for discussion of some of these issues; see also [12]
and references therein).
My research at Caltech revolved mostly around two very different systems of
nonequilibrium nonlinear oscillators. The first system appears as perhaps the simplest nontrivial example: a Duffing oscillator driven by white noise, which models
interaction with the thermal bath, and also by a sinusoidal driving term, which ensures that the oscillator is not allowed to settle into thermal equilibrium. It is the
interplay of nonlinearity and nonequilibrium that gives rise to the surprisingly rich
effects and even allows us to explore certain fundamental questions in Nonequilibrium
Statistical Mechanics.
The nonlinearity is responsible for bistability for some range of parameters: coexistence of two locally attracting steady-states which differ by amplitude and phase
(with respect to the driving phase). Due to the noise, there is a probability of escape from one attracting fixed point1 to another. However, because the oscillator is
driven, this problem differs from the more straightforward case of thermal activation
in a bistable potential well. In our system, the bistability is not a consequence of
a two-well nature of the (time-dependent) potential. This has interesting implications not only on the switching rate between the basins, but also on the probability
distribution of the oscillator in an appropriate phase space as well.
The problem of noise-induced escape from a metastable state has been studied in
various contexts, from nucleation to chemical reactions to switching in nanomagnets.
The escape rate often has the Arrhenius form W ∝ exp(−R/D), where D is the noise
intensity (temperature, in the case of thermal noise). Of central interest both for
theory and experiment is the activation energy R. Starting with the Kramers paper
In this work we will not deal with complicated attracting structures, such as strange attractors,
so the occasional use of the word “attractor” is to be synonymous with attracting fixed point.

[13] much work has been put into calculating R and the prefactor in the escape rate
for various physical systems. From this point of view, it is particularly important to
reveal generic features of the activation energy R, such as scaling behavior with the
system parameters. The onset of the scaling behavior was noticed first for systems
close to thermal equilibrium [14, 15]. Scaling occurs also for systems far from thermal
equilibrium, which generally do not have detailed balance [16, 17, 18, 19, 20]. A
nonequilibrium system may display several types of scaling behaviors and crossovers
between different scaling regimes with varying system parameters. The importance of
studying scaling behaviors and scaling crossovers in a nonequilibrium system comes
from the idea that scaling exponents may be generic and apply to a wide range of
systems [21]. Chapter 3 is devoted specifically to the crossovers between different
scaling regimes in the nonlinear oscillator.
On the experimental front, several types of scaling behavior of the switching activation energy have been observed in recent years for such different experimental
systems as modulated nanoelectromechanical systems (NEMS) [4, 22, 23], and microelectromechanical (MEMS) oscillators [24] (various sources of noise in the nanomechanical systems have been thoroughly described in [25]), Josephson junctions [26],
optically trapped atoms [27], and superconducting resonators [28]. All of these systems are modelled as Duffing oscillators. Not only is this scaling interesting on its
own, but it also underlies various applications, an important example being quantum
measurements [26].
The lack of detailed balance is fundamentally important here: a nonequilibrium
system obeying detailed balance is in many respects similar to equilibrium system.
One remarkable fact about a nonequilibrium system obeying detailed balance has been
proven by Graham and Haken [29], who showed that the probability density function
(PDF) over the phase space coordinates ~r can be expressed in the Boltzmann-like form
exp[−φ(~r)/D], with a globally smooth effective potential φ(~r). The other, closelyrelated manifestation of detailed balance is that the most probable fluctuational trajectory to a given state is related to the fluctuation-free (relaxational) trajectory back
to the stable state by time reversal.

The detailed balance may be fragile with respect to the driving, as indeed happens
in the case of the sinusoidally driven Duffing oscillator. We show that the detailed
balance is broken at any nonzero value of the driving strength. The question remains,
however, whether the smoothness of φ(r) is fragile? Chapter 4 is devoted to this
question and to detailed discussion of detailed balance in our particular system. We
found that in the limit of vanishing noise strength, the tails of φ(~r) become piecewise
differentiable at an infinitesimal driving strength, and we have found the specific
mechanism by which this happens. It turns out that nonsmooth features of φ(~r) arrive
from infinity and move in closer to the attractor of the oscillator as the strength of the
driving is increased. The distance of these singularities from the attractor versus the
driving strength displays power laws with crossovers. Again, this is closely-related
to the pattern of fluctuational trajectories: the most probable paths to neighboring
states are far away from each other, and are not simply-related to the relaxational
paths by time-reversal.
The φ(~r) is the action of the auxiliary problem in classical mechanics that comes
out of the WKB treatment of the Fokker-Planck equation for the PDF. Nondifferentiability of φ(~r) occurs as a result of multivaluedness of the action [30, 31]: only one
of the branches must be chosen in the limit D → 0. Such discontinuities are related
to the objects called caustics - the envelopes of trajectories. Here they are envelopes
of most probable trajectories from an attracting fixed point to a point ~r (in analogy
with Quantum Mechanics they may be called “classical trajectories”). Chapter 4 will
explore the development of nonsmooth features of the PDF as the driving strength is
increased from zero. The formation of regions where the action is multivalued, and
the resultant rich phenomenology exhibited by φ through caustics and switching lines
is ultimately related to the nonintegrability of the classical mechanics that comes out
of the WKB treatment. Understanding the conditions required to maintain integrability in the absence of detailed balance remains at frontier of mathematical physics
[32]. This problem is of fundamental relevance also because the aforementioned mechanism of development of nonsmooth effective potentials is believed to be generic in
nonequilibrium systems described by finite number of degrees of freedom [30, 31].

In summary, even very simple systems lacking detailed balance display a rich phenomenology, and are not well understood; the conditions under which they retain
properties of systems with detailed balance, and hence are qualitatively similar to
equilibrium systems is one appealing line of research. It would also be very interesting to consider spatially extended systems and to explore noise-induced transitions
as a mechanism for pattern selection.

While part I of this dissertation is devoted to nonequilibrium effects in a single
nonlinear oscillator, part II of the thesis is devoted to nonequilibrium collective effects
in a 1-dimensional chain of nonlinear oscillators. One of the most exciting aspects of
modern physics is the investigation of collective phenomena in many-body problems.
Spontaneous synchronization in networks of interacting nonlinear oscillators with different, often randomly distributed frequencies is one of the best understood collective
phenomena in nonequilibrium systems [33]. If the interaction between the oscillators
is weak, then they will evolve in an uncoordinated fashion. In contrast, with proper
connectivity and coupling strength, the population can entrain to the common frequency, a process that can be reflected by an increase of an appropriately defined
order parameter. This collective behavior has been observed in diverse systems including neural networks [34, 35, 36], Josephson junctions [37, 38, 39], lasers [40],
electronic circuits [41], and a range of biological systems, including sleep cycles [42],
animal gaits [43], and biological rhythms in general [44]. Synchronization, as one of
the simplest examples of collective phenomena in nonequilibrium systems, has also
been the focus of theoretical work [45].
Useful insights into the phenomenon of synchronization are obtained by analyzing
the idealized case of all-to-all coupling in the limit of a large number of oscillators. In
one of the most common models each oscillator is described by a single variable, phase,
that advances at a constant rate in the absence of coupling. This model is exactly
soluble using a mean-field approach, and it provides a framework for understanding
synchorinzation [45, 46, 47]. In this model, the synchronized phase in which a finite
fraction of oscillators evolve at a common frequency appears above a critical coupling

strength. Many features of this phase, such as the growth of the fraction of locked
oscillators as a function of the coupling strength, can be calculated.
In many physical systems finite-range interactions depict a more realistic picture.
The simplest of these interactions are nearest-neighbor couplings [48, 49]. In general,
finite-range coupled systems are much harder to analyze, and, as a result, there
are few theoretical tools that produce a quantitative description of the collective
motion. Much of our understanding of systems of coupled oscillators with shortrange interaction rests on numerical solutions of the governing equations. Yet, such
systems are especially important from the fundamental point of view, for it is an
amazing fact of nature that short-range coupling can give rise to collective phenomena
on a much longer, even macroscopic length scale. Understanding of the mechanism of
this process for equilibrium systems is the great achievement of the theory of critical
phenomena [50] and deserves to be better understood for nonequilibrium systems as
well, particularly in non-Hamiltonian models.
In part II we analyze a one-dimensional chain of nearest-neighbor coupled phase
oscillators with quenched disorder in coupling strengths and intrinsic frequencies.
Despite the relative simplicity of the model, it can exhibit complex behavior since
randomness competes with a tendency to establish macroscopic order: Strogatz and
Mirollo [51, 52] nicely showed that a nearest-neighbor chain cannot exhibit extensive synchronization. Nevertheless, the chain develops collective structures, albeit
at finite length scales. These structures are clusters of common average frequencies,
although the instantaneous dynamics of each oscillator in a given cluster can be quite
complicated. The size of these clusters grows with increasing coupling strength.
We propose a method alternative to the direct numerical simulation to attack the
synchronization problem in 1 dimension. The method we develop is a real-space renormalization group (RG) designed to work well on systems with strong randomness,
where both the intrinsic oscillation frequencies and coupling constants are taken from
distributions with long tails. The specific example we investigate contains Lorentzian
distributions. We apply our RG technique to the problem of frequency clusters and
make accurate predictions of various statistical properties of these clusters.

Part I
Aspects of Thermal Activation of
Nonlinear Oscillators out of
Equilibrium

Chapter 2
Driven Nonlinear Oscillator
Interacting with a Heat Bath
This is a background chapter. The concepts developed here will be useful throughout
the subsequent chapters. We will be concerned with a single driven Duffing oscillator
interacting with the heat bath. The starting point for us will be the Langevin equation
description,
q̈ + 2Γq̇ + ω02 q + γq 3 = A cos (ωF t) + f (t),

(2.1)

which is not the most fundamental, but one in which the complicated interaction with
the bath has been reduced to a single additive white noise term f (t)
hf (t)f (t0 )i = 2Bδ(t − t0 ).

(2.2)

The “passive” interaction with the environment is modelled here by a simple damping
and noise term; the steps that lead to this simplified description form is a long and
separate research topic by itself. This reduction has been described for a nonlinear
oscillator in [53] for example. The rest of the terms represent Newton’s second law
for an oscillator in a nonlinear potential and driven via the sinusoidal driving. A wide
variety of systems can be modelled as such. The specific case of a doublyclamped
beam is worked out in Appendix A, where the steps leading from the PDE for the
dynamics of the beam to the Duffing-like equation are shown.
Some of the nontrivial questions that can be asked are: (1) what is the effect of the

noise and what is the probability density function (PDF) ρ(~r) with the presence of
the noise term1 and (2) what are the properties of this nonequilibrium ρ(~r) and how
do they differ from the equilibrium ρ(~r)? The body of physics literature on this topic
is too vast to cite at once; different sources will be cited throughout the discussion. At
present, there remain unanswered questions, some of which will be addressed in part
I of this thesis. From the mathematical point of view, the classic work describing the
effect of noise on dynamical systems is a text by Wentzel and Freidlin, [54]. Another
popular text on stochastic methods in physical sciences is [55].

2.1

From Langevin Equation to Lagrangian Manifolds

To address the big questions posed above we will reduce the problem from that of a
noisy ODE to a calculation of a certain action of a certain “auxiliary” Hamiltonian.
The steps that lead to this auxiliary description are not new, but for clarity and
convenience I choose to briefly go through the steps that lead to this description.
The sources related to various steps would be cited throughout. Readers wishing to
skip over the details, are invited to go directly to the “summary” paragraph with
Eqs. (2.36 - 2.40) and read until the end of Sec. 2.1.1, and then skip further until
section 2.1.3.

2.1.1

Slow Envelope Dynamics

The aim here is to describe the behavior of the Duffing oscillator in a particular
regime when it is weakly driven close to the resonance. In such a regime the oscillator will be nonchaotic, and its dynamic is given by simple harmonic oscillations
modulated by the slowly varying modulation. For a more exhaustive survey of the
Duffing oscillator, and the other regions of its parameter space one is invited to con1

Given a point in phase space at ~r, the quantity ρ(~r)dV is the probability to find a system in
the small neighborhood of ~r of volume dV . In the ensemble picture, this gives the fraction of the
elements of the ensemble that are found in this small region of phase space.

10
sult sources such as [56, 57, 58] and references therein. Thus, our first step is to
derive the amplitude equations describing the slow modulation of simple harmonic
motion due to the influence of the driving, damping and the cubic term in the Duffing
equation. This will be done in the special case when the oscillator is driven close to
linear resonance, i.e., ωF ≈ ω0 . The derivation of the amplitude equations of a nonlinear oscillator is described in many sources, all using somewhat different perturbation
techniques and somewhat different notations; for example Nayfeh [58, 59] does this
with the method of multiple scales and a classic text by Bogoliubov and Mitropolsky
[60] employs their own elegant set of perturbation methods (the Mechanics volume of
Landau and Lifshitz series [62] contains a typically terse, but a physically insightful
discussion of nonlinear oscillations). For completeness I chose to go through the details here employing the averaging perturbation method (this and other perturbation
techniques for dynamical systems can be found in [61] for example). The outcome of
this procedure appears at the top of page 16.
It proves useful to start by nondimensionalizing Eq. (2.1). Let
x =
F =
f˜ =
τ =
λ =
Ω =

γ 1/2
q,
ω0
γ 1/2
A,
ω03
γ 1/2
f,
ω03
ω0 t,
ω0
ωF
ω0

(2.3)
(2.4)
(2.5)
(2.6)
(2.7)
(2.8)

Note that q0 ≡ ω0 /γ 1/2 is the displacement at which the linear and nonlinear terms
are equal, so x is a displacement measured in units of this characteristic length x0 .
Similarly, F0 ≡ ω03 /γ 1/2 is the force needed to be applied to extend q to q0 . In terms

11
of all these parameters, Eq. (2.1) reads
d2 x
1 dx
+ x + x3 = F cos (Ωτ ) + f˜(τ ).
dτ
Q dτ

(2.9)

The Q = ω0 /2Γ is the quality factor which is a dimensionless measure of the damping.
The correlation properties of the noise are modified as well:
hf˜(τ )f˜(τ 0 )i = 2B 0 δ(τ − τ 0 ),

(2.10)

where
B0 =

B.
ω05

(2.11)

We now specialize to the case of large Q, which implies a natural smallness parameter ² ≡ 1/Q, in terms of which we can do perturbation theory. In this regime,
the width of the resonance curve is much narrower then the resonance frequency. Let
us specialize to the physics close to resonance. Set
Ω = 1 + ²σ,

(2.12)

where σ ∼ O(1). We know that near resonance a linear oscillator’s response is
x ∼ F/². Let F = ²p and study what happens as p is varied. We have
F ∼ ²p ,

(2.13)

ẍ, x ∼ ²p−1 ,

(2.14)

x3 ∼ ²3(p−1) ,

(2.15)

²ẋ ∼ ²p .

(2.16)

Clearly we want p > 1 for the linear term to be bigger then the nonlinear one,
otherwise the oscillator is in the strongly nonlinear regime and the linear relationship
x ∼ F/² does not apply in the first place. We would like the force to be not too large
so that the effect of nonlinearity is perturbative, but the linear terms are dominant.

12
When p = 3/2, the cubic, the damping, and the driving terms balance. For p > 3/2
the cubic term becomes much smaller then the damping and the forcing terms, which
are both ²p . So let F ∼ ²3/2 (we set F ≡ F 0 ²3/2 ) and work out the consequences; it
becomes clear that the consequences apply to the cases when F 0 ¿ O(1) which is
equivalent to letting p > 3/2.
Under these conditions the dynamic is approximately that of the lightly damped
harmonic oscillator. In general, it has a modulated envelope dynamics due to the
beating of the natural and driving frequency, but this beating will decay away and
some steady-state amplitude and phase (with respect to the driving) will be reached.
Our objective is to calculate this envelope dynamics and the steady-state amplitude
and phase to which it relaxes. This is where the nonlinearity plays a role; the linear
solutions are known, and the nonlinearity will modify this envelope dynamics and the
steady-state. To obtain an equation for the envelope dynamics we follow a common
method in the theory of nonlinear oscillations. We make a transformation
4/3 (Y sin (Ωτ ) + X cos (Ωτ )) ,
ẋ = Ω²1/2 4/3 (Y cos (Ωτ ) − X sin (Ωτ ))
= ²1/2 (1 + ²σ) 4/3 (Y cos (Ωτ ) − X sin (Ωτ )) .
x = ²1/2

(2.17)

(2.18)

Note that if the oscillations are simple harmonic, the Y and X are constants. If
the motion is not simple harmonic, Y and X are generally functions of time (τ )
which represent precisely the modulated envelope dynamics. They can be found by
substituting Eqs. (2.17 - 2.18) into the equation of motion, Eq. (2.9) and projecting
out the sine or the cosine. We now do just this. First,
1/2
3/4ẍ = Ω²
Ẏ cos (Ωτ ) − Ẋ sin (Ωτ ) − Ω2 ²1/2 (Y sin (Ωτ ) + X cos (Ωτ ))
= ²1/2 (1 + ²σ + ...) Ẏ cos (Ωτ ) − Ẋ sin (Ωτ )
− ²1/2 (1 + 2²σ + ...) (Y sin (Ωτ ) + X cos (Ωτ )) ,

(2.19)

where we used Eq. (2.12). We now substitute Eqs. (2.17 - 2.19) into Eq. (2.9). Some

13
terms cancel.
²1/2 Ẏ cos (Ωτ ) − Ẋ sin (Ωτ ) + ²3/2 σ Ẏ cos (Ωτ ) − Ẋ sin (Ωτ )
− 2²3/2 σ (Y sin (Ωτ ) + X cos (Ωτ )) + ²3/2 (Y cos (Ωτ ) − X sin (Ωτ ))
+ ²3/2 (4/3) Y 3 sin3 (Ωτ ) + 3Y 2 X sin2 (Ωτ ) cos (Ωτ ) + 3Y X 2 sin (Ωτ ) cos2 (Ωτ ) + X 3 cos3 (Ωτ )

31/2 0 3/2
31/2 0 3/2
F ² cos (Ωτ ) +
f²

(2.20)

We pulled out a factor of ²3/2 from the noise function: f˜(τ ) = f 0 (τ )²3/2 . This is done
without loss of generality: it will simply be reflected in the correlation property of
f 0 . We also notice that the O(²1/2 ) terms suggest that Ẋ and Ẏ are independent of
τ . The higher-order corrections more correctly suggest that these terms are simply
small, so on the timescale τ , the functions X and Y vary slowly, reflecting the fact
that the beating near the resonance is slow. It will be useful to define below a new
timescale T 0 = ²τ , which is the characteristic timescale of these functions X and Y .
We multiply the rest by sin (Ωτ ) and integrate over [τ, τ + 2π/Ω] or we multiply by
cos (Ωτ ) and also integrate over [τ, τ +2π/Ω]: this is done while treating the functions
X and Y as constants. We also use the following condition
dY
dX
cos (Ωτ ) = −
sin (Ωτ ),
dτ
dτ

(2.21)

which came from differentiating Eqs. (2.17) and setting it equal to (2.18). After this
procedure, the following equations for Y and X are obtained (we now switch to the
slow timescale T 0 )
1£
dX
−(2σ)Y
(X
dT 0
2·
dY
31/2 0
F + NY ,
−Y + (2σ)X − X(X + Y ) +
dT 0

(2.22)
(2.23)

14
where
Ω τ +2π/Ω 0 0
NX =
f (τ ) sin (Ωτ 0 ) dτ 0 ,
π τ
Ω τ +2π/Ω 0 0
NY =
f (τ ) cos (Ωτ 0 ) dτ 0 .
π τ

(2.24)
(2.25)

We finally make one more round of variable changes. Let η = 1/(2σ), T = T 0 /(2η),
1/2

1/2

Y = P/η 1/2 , X = Q/η 1/2 , F = 3 2 F 0 η 3/2 , Ni = 3 2 Ni η 3/2 and obtain
dQ
= −P − ηQ + P (Q2 + P 2 ) + NQ ,
dT
dP
= Q − ηP − Q(Q2 + P 2 ) + F + NP .
dT

(2.26)
(2.27)

Without noise these would be called the “envelope equations,” or the “amplitude
equations” with which we will work throughout the entire part I of this thesis. While
arriving at them, we went through several rounds of redefining variables, which are
summarized for convenience of the reader in Table 2.1. We track down all of those
renormalizations and express q and q̇ in terms of P and Q and the parameters of the
original Duffing model. The result is

ω0 p
|ωF − ω0 | (P sin (ωF t) + Q cos (ωF t)) ,
r p
dq
23/2
ω0
= 1/2 ωF
|ωF − ω0 | (P cos (ωF t) − Q sin (ωF t)) ,
dt
23/2
q = 1/2

(2.28)
(2.29)

and the two parameters in the amplitude equations are
η =

ωF − ω0

F = A

(2.30)
3γ

32ω03 |ωF − ω0 |

(2.31)

We can also track down all the changes to the noise terms to relate the correlation
functions of the original noise to the correlation functions of the noise in the amplitude
equations. With details, these calculations are somewhat tedious, so they have been

15
relegated to Appendix B. The results can be summarized as follows. Given the
convention
hNi (T1 )Nj (T2 )i = 2Di,j (T1 , T2 ),

(2.32)

where D is the “diffusion matrix,” we have

3γη 2
Di,j (T1 , T2 ) =
B δi,j δ(T2 − T1 )
16ω03 Γ2
= Dδi,j δ(T2 − T1 ),

(2.33)
(2.34)

where
D≡

3γ
16ω03 (ωF − ω0 )2

(2.35)

expressed in terms of the parameters of the natural oscillator defined in the beginning
of this section.

Variable
Position, q
Driving
,A
strength

Table 2.1: Summary of variable changes
Change 1
Change 2
Change 3
X(x, ẋ)
1/2
Y (x, ẋ)
x = γω0 q
Q = η 1/2 X, P = η 1/2 Y
Eqs. (2.17)-(2.18)
1/2

F = γω3 A

1/2

F 0 = ²3/2

F = 3 2 F 0 η 3/2
NQ =

Noise, f
Time, t

1/2
f˜ = γω3 f

τ = ω0 t
2Γ
ω0

Damping, Γ

² ≡

Driving
, ωF
frequency

Ω = ωωF0

f˜

f = ²3/2
T 0 = ²τ

31/2 Ωη 3/2
2π

NP =

31/2 Ωη 3/2
2π

R τ +2π/Ω

f 0 (τ 0 ) sin (Ωτ 0 ) dτ 0

R τ +2π/Ω

f 0 (τ 0 ) cos (Ωτ 0 ) dτ 0
T = T2η

η = 2σ

σ = (Ω−1)

periodicity is removed

16
To sum up, we will now work with the following noisy dynamical system:
Q̇ = KQ (Q, P ) + NQ (T ),

(2.36)

Ṗ = KP (Q, P ) + NP (T ),

(2.37)

where
∂g
− ηQ,
∂P
∂g
KP (Q, P ) = −
− ηP.
∂Q

KQ (Q, P ) =

(2.38)
(2.39)

And g(Q, P ) is a Hamiltonian given by
g = (Q2 + P 2 − 1)2 − FQ,

(2.40)

and η is an effective damping which tends to arrest the amplitude to zero. The Ni s
are stochastic terms which will be modelled as white noise with properties defined in
Eqs. (2.2), (2.32 - 2.35).
In the absence of noise, this dynamical system may have one, two or three fixed
points, depending on the damping η or driving strength F. The amplitude versus
F at fixed η are given by curves such as in Fig. 2.1 that exhibit multivaluedness
for η < ηc . The upper and lower branches correspond to amplitudes of attracting
fixed points (or attractors), while the middle branch is the amplitude of a fixed point
with one stable and one unstable eigenvector (a saddle). Similar curve exists for
phases of fixed points versus F. An initial condition (Q0 , P0 ) will evolve into one
of the attractors. In the (Q, P ) space, there is a separatrix between the basins of
attraction of each of the two attractors. This separatrix is formed from a set of initial
conditions that flows into the saddle via its attractive eigendirection, as demonstrated
in Fig. 2.2. For η < ηc , one of the attracting fixed points coalesces with the saddle via
the saddle-node bifurcation at the endpoints of the trivalued region. These endpoints
are denoted by FB and there are two of these: the low-field FBl and the high-field

17
1.4
1.2

Amplitude

1.0
0.8
0.6
0.4
0.2
0.0
0.0

0.2

0.4

0.6

0.8

Figure 2.1: Amplitude of fixed points versus F for three different values of damping:
η = 0.01ηc , η = 0.5ηc , and η = ηc = √13 , left to right.
FBh . As F is slowly varied,2 the system will adiabatically follow one of the attracting
fixed points until this fixed point disappears by colliding with an unstable fixed point
in the saddle-node bifurcation. Upon further varying F the system will settle into
the only other remaining fixed point. Note that reversing the change in F will result
in the hysteretic behavior. Incidentally, there are interesting effects that arise when
F is varied fast. I explored this recently in [63].
The entire η-F parameter space is represented in Fig. 2.3. The closed triangularshaped domain is the region in which there are three solutions. If the parameter range
is in the bistable regime the noise can overthrow the system from one basin to another.
It is critically important to note that at any finite η, bistability will exist only at a
nonzero F. Therefore the problem of transition rates between two stable attractors
is a nonequilibrium problem! Noise-induced transitions from low amplitude to high
amplitude branch of 2.1 will be called “down → up transitions,” and those from high
amplitude to low amplitude branch, will be called “up → down transitions.” In the

Slow variation means that the rate at which the position of the quasi-fixed point changes as F
is varied is much smaller then the rate of relaxation to the attractor determined by its eigenvalues
at that instantaneous value of F and η.

-2

-1

-2

Figure 2.2: The dynamics in the P -Q space for a particular values of parameters
η = 0.3σc ≈ 0.173 and F = 0.3. The fixed points are shown as black solid dots. The
left-most one is the saddle, which lies on the separatrix between two basins. The other
curves shown are trajectory from the saddle to each of the attractors (blue spiraling
curves) and some other typical trajectory (dashed curve).
limit of very small noise, when D is the smallest parameter in the system, there is a
characteristic time for these switching events. In this limit, the switchings are rare
events because the system has to fight the much stronger tendency of an oscillator to
return to one of the fixed points. The escape process is cumulative and can be thought
~ in Eqs. (2.36)-(2.37). A
of as a diffusion process in the deterministic drift field K
mathematical description of this process is captured by the Fokker-Planck equation
which we describe next. Intuitively, a system placed at one of the attractors at time
t0 will wonder around this attractor and occasionally find itself far enough away from
it to reach and overstep the separatrix. Due to the strong deterministic current, the
most likely subsequent scenario will be the quick relaxation of the system toward the

19
0.6
0.5

0.4
0.3
0.2
0.1
0.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

Figure 2.3: The region of bistability in the η-F parameter space is the closed region
between
the FBh (top)
q and
q the cusp point
´ FB (bottom) curves. These curves meet at

. The width of the hysteresis at zero η is 27
. Outside of
at ηc = √13 , Fc = 27
the bistability region there is always only one fixed point which is an attractor.

attractor on this new side of the separatrix, completing the switching event.

2.1.2

Fokker-Planck equation and the Path Integral solution

The probability density function (PDF) ρ(Q, P ) of reaching (Q, P ) from some initial
(Q0 , P0 ) satisfies the Fokker-Planck equation (FPE), see [30] for example:
∂ρ
= −∇ · Kρ − D∇ρ
∂t
~ · ~j.
= −∇

(2.41)
(2.42)

This can be obtained by standard derivation methods (see [2] or [64] for example),
but it has a simple intuitive interpretation of diffusion in a drift. The ~j is simply a
probability current, which has a deterministic component and a diffusive component,
much like ink dropped in a current will be carried away by the current, and it will also
spread out diffusively. This FPE can be rewritten to look like a Liouville equation

20
with a scattering term on the right-hand side:
∂ρ
− {g, ρ} = ∇ · η~rρ + D∇ρ
∂T

(2.43)

The solution can be written as a functional integral
ρ(~r, t|~r0 ) =

Dµ [~r] exp −

Z ~r=~r(t)

L ~r(τ ), ~ṙ(τ ) dτ ,

(2.44)

r0 =~
r(0)

over all paths [30, 54, 65] where ~r represents a position in the (Q, P )-space, Dµ[~
ρ] is
the functional measure and the Lagrangian is
´2 ³
´2 ¸
1 ³
L=
Q̇ − KQ + Ṗ − KP
+ O(D).

(2.45)

The O(D), and an even smaller O D3/2 terms will not be considered: we are currently interested in the dominant contribution in the limit of vanishing noise strength.
The leading part of L that we consider has a very simple geometric interpretation.
The probability density of a given trajectory can be computed by first breaking up
this trajectory into small pieces, each piece corresponds to a constant time interval
~ is constant. Suppose the prob∆t. Consider the kthe piece: locally the drift field K
ability of being found at the beginning of this kth piece at time tk is known (i.e.,
the PDF is a delta function of unit strength). Over the course of the subsequent
∆t, the probability density will spread out and it will also drift. We can then use
an elementary solution to the diffusion problem to calculate the probability of being
found at the end of the kth piece after time ∆t, i.e.,
³
´2 
~ṙ∆t − K∆t
ρk ∼ exp 
,
4D∆t

(2.46)

~ are evaluated at position of the beginning of kth piece, ~rk . Finally,
where ~ṙ and K
the probability of reaching ~r from ~r0 is the product Πk ρk : the product of exponentials
gives an exponential of a sum. We also see that since the difference between ~ṙ and

21
~ is just the noise, the leading-order part of the action represents the integral over
¯2
R ¯¯
~ (T )¯¯ dT , which can
the noise realization which allows reaching ~r and minimizes ¯N
be interpreted as the “noise energy,” or the work the system needs to do against the
stochastic force. Thus, the sum over trajectories can also be represented as the sum
over noise realizations.
A rigorous and general derivation by Graham can be found in [66] (for a concise
summary, see Section III of Graham’s 1985 paper [30]), which also includes the omit¡
ted O(D) and O D3/2 terms, as well as a generalization to the case of a nonunit
diffusion matrix; it also mentions prior works of Onsager and Machlup who first proposed a path integral solution to a PDF, albeit restricted to the Gaussian process. A
complimentary set of sources are those by Dykman [67] and [68] as well as in Dykman
and Smelyanskiy [69], all of which have a more physically intuitive discussions then
those of Graham, but they gloss over some mathematical details which can be found
in the works by Phythian [70] and Jouvet and Phythian [71]. The classics on path
integrals by Feynman and Hibbs [72] and Schulman [73] also serve as useful references
in setting up the path integral.
One can reexpress all paths as variations around a path that minimizes the exponent. We will call this the “optimal path3 .” Then
Z ~r=~r(0) ³
L ~r(τ ), ~ṙ(τ ) dτ
ρ(~r, t|~r0 ) = ρ0 (~r, ~r0 , t) exp − min
r0 =~
r(−t)
S(~r)
= ρ (~r, ~r0 , t) exp −

(2.47)

where the prefactor ρ0 comes from performing path integrals over variations around
the optimal path. This prefactor is generally nonexponential: we can motivate it by
approximating the functional integral by a finite number of regular integrals, each
of which is over a variation and in general converges to some finite quantity. There
are notable exceptions, particularly at the cusp of the pair of caustics (Chapter 4) as
recently pointed out by Maier and Stein in [19]. We will only be concerned with the

In analogy to the path integral formulation of Quantum Mechanics, this can be called the
classical path.

22
exponential contribution, because it will dominate transition rates in the limit of weak
noise. The problem of prefactors for continuous systems without detailed balance
remains open. The formulation of the type in Eq. (2.47), where S(~r) is expressed
as the sum over trajectories ~r(t) is known as the Wentzell-Freidlin formulation. A
different formulation, where the functional integration is also performed over the
noise trajectories and over a Lagrange multiplier that relates noise trajectories and
trajectories r(t) can be found in [67]-[72]. Clearly, the Wentzell-Freidlin functional can
be obtained when noise and the lagrange multiplier trajectories are either integrated
out or a saddle-point type method is used.
We are concerned with the distributions and escape rates in the limit of vanishing
noise, D → 0. Clearly, the escape rates W from one basin to the other approach ∞ in
this limit.4 In an ensemble of experiments where all systems are initially positioned
at some ~r0 , the distribution ρ(~r, t|~r0 ) quickly establishes a quasistationary form on
the time scale τrel ¿ t ¿ W −1 , where τrel is the relaxational time scale dictated
by the noise-free dynamics. In the limit D → 0 we can replace this time range
by τrel ¿ t < ∞. Then, we can say that in the limit D → 0 the quasistationary
distribution ρ(~r) is given by limt→∞ ρ(~r, t|~r0 ) (note that D = 0 and the limit D → 0
are not the same). This limit is indeed independent of ~r0 and hence represents the
steady-state probability distribution [30]. Mathematically, −D ln [limt→∞ ρ(~r, t|~r0 )]
is the action of a trajectory from ~r0 to ~r that takes an infinitely long time. With
no finite time constraint to reaching ~r, the most likely trajectory from ~r0 to ~r first
~ (i.e. the noise-free path governed by
follows a relaxational path to the attractor of K
~ this costs zero action, and takes an infinitely long time, so all the memory
~ṙ = K);
of ~r0 is lost, because all such trajectories to reach ~r will pass through the attracting
~ corresponds to zero-momentum
fixed point. Notice also from Eq. (2.45) that ~ṙ = K

Note that there are two contexts in which the word “escape” appears. (1) escape from an attractor to a point ~r: this most likely happens via an optimal (or classical) trajectories and moreover, the
energy of this trajectory is 0 in the case of stationary states. These may also be called fluctuational
trajectories, as opposed to the relaxational ones. (2) escape from one attractor to another, as used
in this sentence and in Sec. 2.1.4: this most likely happens through a special optimal trajectory
which crosses the separatrix through the saddle (see below). Reaching a saddle normally involves a
large deviation from the attractor: in our work the noise is assumed to be the smallest parameter
in the theory, so this is always a rare event.

23
“auxiliary” dynamics of the Lagrangian L. The Hamiltonian associated with this
Lagrangian L is

H = p2Q + p2P + pQ KQ + pP KP .

(2.48)

This H is often called the “auxiliary” Hamiltonian, to avoid confusion with g. The
points (QF P , PF P , pQ = 0, pP = 0) are fixed points of the auxiliary dynamical system
formed from H. Therefore, the remaining part of the path, from the attractor to ~r,
also takes an infinite time because it escapes the fixed point of the auxiliary dynamics;
however, unlike the first part of the path, it will cost a finite action.
In reality, the quantity −D ln [limt→∞ ρ(~r, t|~r0 )] is the leading-order term in the
expansion of the action in powers of D and the resulting ρ represents stationary
distribution. The corrections in powers of D would take into account the fact that
W −1 is finite. In short, the leading-order exponential part (the stationary part) of
the quasistationary distribution is given by

ρ(~r) ∼ exp − min

2.1.3

Z ~r=~r(0)

L ~r(τ ), ~ṙ(τ ) dτ .

(2.49)

ra (−∞)

Auxiliary Classical Mechanics

To calculate ρ(~r) up to this accuracy in D, we need in principle to solve EulerLagrange equations to find a path that leads from the attractor to point ~r = (Q, P )
and evaluate the action S along this path. The action S(~r) serves as a generalized
potential. Note that there are infinitely many optimal (or classical) trajectories from
an attractor to the point ~r; they all differ by their initial momenta, but only those
trajectories that have zero energy arrive to and emanate from fixed points (in 4D)
and thus take an infinite time. In fact, another way to see that only zero energy (and
hence infinite time) trajectories are related to stationary ρ(Q, P ) is to substitute
ρ(Q, P ) = ρ0 (Q, P ) exp (−S(Q, P )/D) into the stationary FPE [74], and obtain

∂S
∂Q

¶2

∂S
∂P

¶2

∂S
∂S
KQ +
KP = 0.
∂Q
∂P

(2.50)

24
∂S ∂S
This is just a Hamilton-Jacobi equation, H(Q, P, ∂Q
, ∂P ) = 0, and it describes dy-

~ F P , PF P ) = 0, the momenta at all
namics on the H = 0 manifold. Now, since K(Q
the FPs must also be zero.
The Hamiltonian EOM that follow from H form an auxililary 4D dyanmical system
with two momenta, pQ and pP and two spatial variables Q and P . This system is
∂g
= 2pQ − ηQ − P + (Q2 + P 2 ), P
(2.51)
∂P
∂g
= 2pP − ηP + Q − (Q2 + P 2 )Q + F,
(2.52)
Ṗ = 2pP − ηy −
∂Q
∂ 2g
∂ 2g
ṗQ = ηpQ − pQ
= (η − 2QP )pQ + (3Q2 + P 2 − 1)pP , (2.53)
+ pP
∂Q∂P
∂Q
∂ g
∂ 2g
ṗP = ηpP − pQ 2 + pP
= −(Q2 + 3P 2 − 1)pQ + (η + 2QP )pP .(2.54)
∂P
∂Q∂P
Q̇ = 2pQ − ηQ +

We notice immediately that the deterministic 2D dynamics takes place on the p~ = 0
~ = noise, so in the absence of noise, the only possible
plane. In fact, 2~p = 21 (~ẋ − K)
dynamics lie on the p~ = 0 plane. Any fluctuational (as opposed to relaxational) dynamics will be lifted out of this plane. As already mentioned, the fixed points (FP)
of this dynamics are (QF P , PF P , 0, 0) where (QF P , PF P ) are the FP of the original
2D dynamics. It can also be easily shown that eigenvalues of each of these auxiliary
FP are (λ1 , λ2 , −λ1 , −λ2 ) where λ1 and λ2 are eigenvalues of the FP of the original
dynamics. The p~ = 0 plane forms the stable manifold of each of the fixed points. A
trajectory that escapes from an attracting fixed point (now using the first meaning
word escape, see footnote 4) will lie on a certain manifold (surface) in the auxiliary
(Q, P, pQ , pP ) space. In the vicinity of the attracting fixed point, this manifold is
tangent to the plane spanned by two unstable eigenvectors of this FP. Hence, the
manifold on which the fluctuational trajectories lie is conventionally called the unstable manifold. The stable and unstable manifolds are known as Lagrangian Manifolds
[30, 74, 75, 76]. Calculation of trajectories that escape the attractor will involve
placing initial conditions on this unstable manifold close to the attracting FP. This
procedure is described in detail in the following Chapter. Some analytical aspects of
the Classical Mechanics associated with H in chapter 4.

2.1.4

Escape Rates

To calculate the escape rate [65, 67] from one attracting fixed point to the other (here
I am making use of the second meaning of the word escape), in principle, one must
compute the line integral of ~j(l) · n̂(l) over the separatrix, where ~j is the probability
current, n̂ is the unit normal to the separatrix and l is the arc length along the
separatrix. The ~j must come from a nonstationary distribution, otherwise the net flow
from the basin would be zero in view of divergence theorem (see Eq. (2.42) with zero
left hand side). For example, it may come from the slowly-decaying quasistationary
distribution. We have seen that the leading order term in D of such a distribution
is given by the stationary distribution in Eq. (2.49). Moreover, it may be seen from
the saddle-point method [30], for example, that the line integral over the separatrix
will pick up an exponential contribution given by the minimum value of the integrand
along the separatrix. In fact, it was shown by Dykman and Kryvoglaz [65] that of all
the classical trajectories emanating from the attractor and leading to the separatrix
around its basin of attraction, the path that leads to the saddle point, if it lies on
the separatrix (as it does in our system), will have the minimum action. We call this
action S ∗ . Furthermore, in performing the line integral over the separatrix, one can
factor out exp (−S ∗ /D) and perform the integration, which will give some number.
In other words, up to the exponential factor, the escape rate W will simply be given
by
W ∼ exp (−S ∗ /D),

(2.55)

which also happens to be the exponential factor in ρ(~rs ). The preexponential factors
(of either the distribution ρ(~r) or escape rate W ) are not easy to obtain but they are
also the least important in the limit of low D!
In summary, in the limit of weak noise, the rate of escape W from one attracting
fixed point to another is given by W ≈ exp (−R/D) where R = S ∗ is the action along
the trajectory which connects the attracting fixed point to the saddle and satisfies
H = 0.

2.2

Analytical Limits

Both limits discussed below are physically analogous to the two types of limits considered by Kramers [13]: an overdamped limit and a low-damping limit. One major
difference is that our Hamiltonian g(Q, P ) does not have an explicit separation of
kinetic and potential energies. The results in Chapter 3 will be compared with the
analytical limits derived in this section.

2.2.1

Close to Bifurcations

At either bifurcation point (see Fig. 2.1) which we called FBl or FBh , the saddle and one
~ merge together at (QB , PB ) in a saddle-node bifurcation.
of the attracting nodes of K
As F approaches the bifurcating value, FB , a saddle and a node approach each other,
the attractive eigenvector of the saddle and one of the attracting eigenvectors of the
node align parallel to each other, and the corresponding attracting eigenvalues of both
FP approach the same negative value. Also, the repulsive eigenvector of the saddle
and the other attracting eigenvector of the node become equal and the corresponding
eigenvalues both tend to zero. This process is depicted in the cartoon in Fig. 2.4.

Figure 2.4: Cartoon of the K-flow
near a saddle-node bifurcation.
Thus, close to a bifurcation, there develops a soft mode: a narrow region connecting
the saddle with an attractor along which the K-motion
is slow. Due to this slowing
down, the soft mode forms a path of least resistance along which a large noise-

27
induced fluctuation away from the attractor is most likely to take place. To analyze
~ in terms of parameter δ = |F − FB | and variables (f =
this system we rewrite K
Q − QB , s = P − PB ); it turns out to be unnecessary to make a linear transformation
~ at (QB , PB , FB ). In terms of these new coordinates, K
to eigencoordinates of K
transforms to the following form:
f˙ = −2η(f − aB s) + PB f 2 + 2QB sf + 3PB s2 + f 2 s + s3 + NQ

(2.56)

ṡ = ±δ − 3QB f 2 − 2PB f s − QB s2 − f 3 − s2 f + NP

(2.57)

Here the + sign applies for the low field (high amplitude) bifurcation and
aB = η −1 (1 ± 2 1 − 3η 2 )/3.

(2.58)

Notice that the position of the fixed points sets a characteristic scaling of both s and
f variables to be ∝ δ 1/2 . From this we see that s is a slow variable while f is fast on its timescale, s appears approximately frozen. To lowest order in δ, we treat s as
completely frozen, in which case f relaxes to aB s plus corrections of order δ. In such
adiabatic limit, the dynamics of s is then given by
ṡ ≈ ±δ − b(η)s2 + NP ,

(2.59)

b(η) = 3QB a2B + 2PB aB + QB ,

(2.60)

where

(b > 0 for the low-field bifurcation and b < 0 for the high-field bifurcation). Notice
that the characteristic timescale for the variable s, characterized by the eigenvalues
around the fixed point, scales as δ 1/2 , while the relaxational timescale for the variable
f scales as 2η which is not assumed to be particularly small, so indeed s is a slow
variable. The equation describing the dynamics of the slow variable, Eq. (2.59), can

28
be seen as an equation for an overdamped particle in a cubic potential U given by
U = ∓δs + s3 ,
while the distance between the saddle and the attractor is

(2.61)
δ/|b|; thus the potential

barrier for the slow coordinate is
∆U =

4δ 3/2
3|b|1/2 (η)

(2.62)

This is the well-known ξ = 3/2 activation scaling. A 1-dimensional Fokker-Planck
equation for s can be written, a stationary solution for ρ(s, η̃, η) obtained, and
the escape rate calculated in the standard way. This quantity will have the form
W = W 0 (η) exp (−∆U/D). This type of analysis for a multidimensional nonlinear
dynamical system was first performed by Dykman and Krivoglaz in 1980 [16]. These
authors also considered the escape problem near the double-bifurcation point in which
the top and bottom branches in Fig. 2.3 meet. Such an adiabatic approach is valid
only sufficiently close to the bifurcation point. When η ¿ 1, “sufficiently close”
means that δ ¿ η 3 . In this regime,
 √
 4 2 η 1/2 δ 3/2 low field bifurcation,
R≈
 4 ηδ 3/2
high field bifurcation.

(2.63)

31/4

A more systematic treatment of this regime will allow us to extract corrections to the
adiabatic limit and derive the δ ¿ η 3 criterion will be considered in Section 3.2.1.2.

2.2.2

Low Effective Damping Limit

~ is approximately Hamiltonian:
When the effective friction η is low, the dynamics of K
the timescale for the decay of energy is much larger then the timescale to make one
cycle on the contour of approximately fixed energy. This separation of timescales
breaks down very close to a bifurcation point (refer back to Figs. 2.1 and 2.3), so
as long as η 6= 0, there will always be a small fraction of the hysteresis displaying

29
a ξ = 3/2 scaling. Aside from this very narrow region, the separation of timescales
allows one to turn a 2-variable FPE into a 1-variable FPE for the diffusion of energy,
by averaging out the fast part of the dynamics. The noise strength D is proportional
to η 2 (see Eq. 2.33), so the FPE in the stationary regime can be rewritten as
~ · ~rρ + D0 η ∇ρ
{g, ρ} = −η ∇

(2.64)

Note that when η = 0, any function of the form ρ(g(Q, P )) will be a solution to this
FPE. Therefore, we can expand
ρ(Q, P ) = ρ(0) (g(Q, P )) + ηρ(1) (Q, P ) + O(η 2 ),

(2.65)

and integrate over a contour of constant g(Q, P ). We will sometimes refer to such
theory as “quasi-Hamiltonian” theory. After some algebra, this leads to:
dE

dρ(0)
B(E)ρ (E) + D ηD(E)
dE
(0)

= 0,

(2.66)

where B(E) and D(E) are given by
Z Z
B(E) =
Z Z

dQ dP,

(2.67)

∇2 g(Q, P ) dQ dP.

(2.68)

g(Q,P )=E

D(E) =
g(Q,P )=E

This approach has been followed by Dmitriev and D’yakonov in [77] An alternative
approach followed by [21] is to calculate an average rate of energy decay due to friction
and an effective diffusion coefficient for the energy drift due to noise. These would
turn out to be precisely the B(E) and D(E) respectively, hence Eq. (2.66) is the FPE
for the energy. Yet an entirely different approach has been undertaken in the paper by
Dykman [65]. Returning back from D0 to D, the activation rate is W ∝ exp (−R/D)
where

Z Es
R = −η
Ea

B(E 0 )
dE 0 .
D(E )

(2.69)

30
Here Ea is the energy of attractor (it is strictly speaking a center at zero η but
becomes an attracting fixed point for finite η), and Es is the energy of the saddle
(it remains a saddle even at zero η). To compute Rup→down the relevant “attractor”
is the one with larger amplitude and to compute Rdown→up the relevant “attractor”
is the one with smaller amplitude (see Fig. 2.1).5 We plot both of these R/η in
Fig. 2.5. It is worth mentioning that in steady-state the ratio of rates gives the ratio

RΗ
0.4
0.3
0.2
0.1
0.2

0.4

0.6

0.8

1.0

4 27

Figure 2.5: Activation barriers in the limit of η → 0 given by Eq. (2.69). The curve
that increases from left to right corresponds to Rup→down and the curve that decreases
from left to right corresponds to Rdown→up .
of occupation probabilities of high or low amplitude states. Thus the F at which
Rdown→up = Rup→down is the point at which the most probable states change. Both
are experimentally relevant quantities and have in fact been measured; see recent
NEMS experiment of Aldrige and Cleland [22] for example.
Close to bifurcation points

R≈

 2ηδ

31/4

low field bifurcation [65],
ηδ 3/2 high field bifurcation [77],

(2.70)

For these down → up transitions there are two disjoint contours of g at a given energy E. The
smaller contour immediately around the low-amplitude center needs to be chosen in the calculation
of B(E) and D(E) in Eqs. (2.67) and (2.68).

31
Of course these equations in particular, and the whole low-damping method (or quasiHamiltonian method) of this subsection breaks down sufficiently close to bifurcation
points when η is finite (because the trajectories are no longer approximately Hamiltonian trajectories), and must be replaced by the strong-damping method (or quasi1-dimensional method) of the previous subsection. As η → 0 the interval of δ over
which the latter is applicable shrinks (as a reminder, δ is the deviation of F from
the value of F at which there is a saddle-node bifurcation point). There must be a
crossover from one regime to the other at some δ. Indeed, comparison of Eqn. (2.70)
with Eqn. (2.63) reveals that there must be a crossover from a ξ = 1 regime to ξ = 3/2
regime as F approaches FBl for up → down transitions near the low-field bifurcation
point. This is the part of my research which is presented in the next chapter.

2.3

Power Law Scalings of Escape Barriers

The aim of this section is to discuss conditions when a generic system is expected
to exhibit the scaling exponents 3/2 and 1. The discussion may appear somewhat
detailed, but for the most part, we summarize recent work of Dykman, Schwartz and
Shapiro [21] in which the scaling exponent at nonzero η was related to the geometry
of homoclinic orbits at zero η.
We just saw that the two analytical limits, the low-η limit and the overdamped
limit close to bifurcation points exhibit their own scaling laws. In general, the inquiry
into scaling exponents is an interesting one because these exponents may be universal.
For example, the 3/2 scaling is ubiquitous, and it also appears in the original work
of Kramers [13]. This is a result of a universal structure of a dynamical system near
the saddle-node bifurcation point. When η 6= 0, sufficiently close to the saddle-node
bifurcation, the dynamics near the attractor and the saddle is generically overdamped
[78, 79]. This leads to the soft mode 1-dimensional geometric picture offered in Section
2.2.1, which leads to the cubic potential, and finally the exponent 3/2. Hence, ξ = 3/2
is a genetic property in a system with finite damping and close enough to the saddlenode bifurcation.

32
Further from the bifurcation, the overdamped soft mode does not exist. We have
seen that for sufficiently low η, the overdamped ξ = 3/2 crosses over to ξ = 1 (close to
FBl ) or remains ξ = 3/2 (close to FBh ). As recently proved by Dykman, Schwartz and
Shapiro [21] (DSS), these are in fact also generic scaling laws in the underdamped
vicinity of a saddle-node bifurcation. DSS considered two scenarios which may happen
when η equals zero. As δ → 0, the center and the saddle may merge, or they may
remain separate. DSS named the first type of Hamiltonian bifurcation local and the
latter nonlocal. The two situations dictate the form of the Hamiltonian g.
Let us switch to coordinates (q, p) such that when η = 0, both the center and the
saddle fixed points lie on the p = 0 axis. In the case of a local Hamiltonian bifurcation,
in the vicinity of (q = 0, p = 0, δ = 0), g is given by gL (q, p, δ) = 21 p2 + Ucub (q) where
Ucub (q) = − 13 q 3 + η̃q. The center lies inside a small homoclinic loop at the energy of
the saddle. DSS proved that given this form of the Hamiltonian g, addition of noise
and damping leads to R(δ) = A(η)δ 3/2 (in the underdamped regime). But we know
that as as the system approaches the bifurcation even closer, it becomes overdamped,
a soft mode sets up and the method of analysis presented in Section 2.2.1 leads to
R(δ) = B(η)δ 3/2 . Therefore, assuming smooth dependence upon δ, it must be that
A(η) = B(η). In other words, the ξ = 3/2 power law in this situation holds beyond
where it was expected to hold due to a soft mode theory! In the Duffing system,
a nonlocal bifurcation takes place at FBh . In the numerical experiments described
below, we found that the ξ = 3/2 power law holds near FBh over a finite range of δ
even as η → 0. This ξ = 3/2 law at zero η has been predicted by [77], but the new
work of DSS gives it a new explanation. We elaborate upon this in Section III.B.
In the case of a nonlocal Hamiltonian bifurcation, g is given by gN L (q, p, δ) =
1 2
p + δU (q), where U (q) is independent of δ and has a local minimum and maximum.

The homoclinic loop of such systems has the property that only its size in the pdirection shrinks to zero as δ → 0, while the distance between the center and the
saddle inside the loop remains finite. Introduction of damping into these types of
systems will not change the separation between the saddle and the attractor until
the system has come close enough to the bifurcation: δ ∼ η 2 . DSS proved that given

33
this form of the Hamiltonian g, addition of noise and damping leads to R(δ) ∝ δ
(in the underdamped regime). Therefore, at finite η we expect a ξ = 1 scaling
moderately close to the bifurcation, but further away then the region of overdamped
ξ = 3/2 scaling. In the Duffing system, a nonlocal bifurcation takes place at FBl . It is
important to stress that all saddle-node bifurcations at finite η are local in the sense
that the attractor and the saddle merge as δ → 0. The terms local or nonlocal only
make sense at zero η. We must also mention that a nonlocal bifurcations may take
place in many forms. For example, in a Duffing oscillator, the homoclinic loop has a
horseshoe shape at η = 0 which becomes thinner while remaining extended as δ → 0.
DSS found a change of variables which maps from g to gN L .

Chapter 3
Scaling Crossovers
In this chapter we study properties of crossovers between different scaling regimes.
Each scaling regime and the corresponding scaling exponent may have a specific physical explanation; we have seen two examples: the overdamped soft mode close to the
bifurcation or the separation of the damped and the hamiltonian timescales further
away from the bifurcation. In such cases an analytical approximation may be possible,
as in fact it was in the aforementioned cases. These analytical approximations often
take the form of a perturbation series. The regime of crossovers between different
scaling exponents is also interesting, however, it is less likely to afford an analytical
approximation because this is precisely where the perturbation theory applicable in
each regime must fail. Perturbation series may still be useful in estimating the characteristic beginning point or the characteristic end point of a crossover region. To
do this one may develop the perturbation series in each scaling regime to one higher
order and look at the values of δ where the correction to the leading order becomes
comparable in magnitude to the leading order itself. We will in fact do this below,
but clearly, this approach is not a rigorous, but one designed to make an estimation.
Thus in understanding crossovers, the power of numerics is very useful and is an
appropriate first step in trying to understand crossovers between scaling regimes.
We would like map out the (η, F) parameter space (cf. Fig. 2.3) for regions in
which each scaling regime holds and for locations of crossover regions between these
scaling regimes. This will tell us the characteristic location of a crossover and the
characteristic width of the crossover region. In what follows we first outline the

35
numerical method used, then present the results and discuss the findings.1

3.1

Numerical Method

Numerical calculations of escape rate are based on direct evaluation of S ∗ , the quantity
defined in Sec. 2.1.3. The algorithm for doing so is as follows. (1) Choose parameters
η and F. (2) Calculate repulsive eigenvectors (i.e., those whose eigenvalues have
positive real part) of the FP (Qattr , Pattr , pQ = 0, pP = 0) of the auxiliary system.
These eigenvectors span the tangent plane to the unstable manifold of this attractor.
(3) Position a locus of initial conditions on a small circle that lies on this tangent plane
and centered around this FP and evolve each of these initial conditions according to
the Hamiltonian equations with Hamiltonian H from Eqn. (2.48). The angle around
this circle, ϕ, serves as a parameter that enumerates a trajectory. 2 (4) From this
set of trajectories, find the one which leads into the saddle, with a special trajectory
parameter ϕ∗ . Because the method involves “shooting” many trajectories to see
which one hits the saddle, this method is sometimes called “the shooting method”.
(5) Evaluate the action along this trajectory. For each η, this procedure would be
repeated for many values of F between FBl and FBh .
In practice, step (4) may require several rounds of bracketing. In each round, we
shoot a given number of trajectories in some range of ϕ, starting from the 2π range
in the first round. For each trajectory we calculate the arc length on the separatrix
between the saddle and the point of intersection of that trajectory with the separatrix.
We identify a ϕ for which the trajectory came closest to the saddle and then evolve
another set of trajectories in a small range of trajectory parameters around this ϕ.
After more and more rounds of such procedure, we obtain a trajectory that hits
the separatrix closer and closer to the saddle. In general, finding a trajectory that
connects an attractor with the saddle is difficult because it lies on the intersection of

A preprint of our paper related to this chapter can be found in [80].
It turns out to be not the ideal way to parametrize trajectories, since the local pattern of spiraling
trajectories may be elliptical, so different values of ϕ may actually lie on the same trajectory. This
made this method a little inconvenient, but it could still be used. A different technique is to sample
initial conditions situated on a straight line in between two windings of a spiral.

36
the unstable manifold of the attractor and the stable manifold of the saddle in the
4D auxiliary space [74]. 3 A small deviation of ϕ on opposite side of ϕ∗ will produce
trajectories that lie on opposite sides of this intersection and will diverge from each
other. This effect is most pronounced closer to bifurcations because the motion along
the optimal trajectory becomes much slower then the diverging motion away from it such is the structure of eigenvalues of the attractor-saddle pair close to the bifurcation.
This polarity takes place only close to bifurcation points, but it persists farther away
from bifurcation points when η is larger. Therefore, in the small-η regime, finding a
trajectory between an attractor and the saddle is difficult (i.e., requires many rounds
of bracketing) only very close to bifurcation points, but in the large-η regime, finding
such a trajectory is difficult over a larger part of the hysteresis.
We will now present results of the computation of S ∗ , scaling laws, and crossovers.
It is sometimes easier to present results in terms of
x≡

FBh − FBl

(3.1)

the “reduced” δ. For each friction η and each type of transition (there are two: up →
down and down → up), we scanned the hysteretic region at multiple values of x, and
for each x we followed the shooting procedure, and where necessary, the bracketing
method. We first summarize our findings and then discuss the details for each type
of transition separately.

The question of existence of this intersection, to our knowledge, is not a fully understood one.

3.2

Results

We first summarize the main results and then discuss them in detail. For the up →
down transitions, we identified three scaling regimes of S ∗ versus F at fixed η: the
power law with exponent 3/2 close to the bifurcation, followed by a regime with the
scaling exponent 1 further from the bifurcation, followed by a nonpower law regime.
The regions of the (η, F) parameter space for each regime are mapped out, showing
the disappearance of the 3/2 regime as η → 0. We point out the consistency of this
behavior with generic predictions of DSS. We also find that for finite η the exponent
3/2 holds over a much wider range of parameters than where it would be expected to
hold based on the overdamped 1-dimensional reduction close to bifurcation points.
For the down → up transitions, we identified a new intermediate scaling with
the exponent of ≈ 1.3, which follows the power law with exponent 3/2 close to the
bifurcation; the nonpower law regime again takes place furthest from the bifurcation.
The regions of the (η, F) parameter space for each regime are mapped out, showing
the persistence of the 3/2 regime as η → 0 (again underscoring the nonoverdamped
nature of this exponent). These observations are also consistent with the generic
predictions of DSS.

3.2.1

Up → Down Transitions

3.2.1.1

Numerical Results

A plot of S ∗ /η versus x, Figs. 3.1, 3.2 reveals several features. At low η two regimes
clearly stand out: one with ξ ≈ 3/2 (regime A), and one with ξ ≈ 1 (regime B).
Further from the bifurcation, the ξ = 1 scaling law breaks down and S ∗ /η versus x
displays a nonpower-law behavior. For various values of η, the data in the regime
A fits well to a power law with ξ which is up to 3% below the theoretical value of
3/2. In general, precise extraction of the exponent 3/2 requires approaching very
close to the bifurcation. It appears that the discrepancy was highest in those cases
when the bifurcation was approached less closely (relative to the crossover point) or

Figure 3.1: Main plots: S ∗ (x)/η (dots) and fits (lines) for up → down transitions.
The fit to the ξ = 3/2 regime was made by analyzing ln (S ∗ /η) versus ln x and a
fit to the ξ = 1 regime was made by analyzing S ∗ /η versus x. Inserts: S ∗ (x)/η
(dots) and theory (lines) based on Eq. (2.63) for the ξ = 3/2 regime and on the
quasi-Hamiltonian theory, Eq. (2.69).

Figure 3.2: S ∗ (x)/η on a linear scale. The thick black line represents the zero damping
theory, Eq. (2.69), the same as one of the curves in Fig. 2.5. The other curves
represent numerical calculations of S ∗ : squares for η = 0.03 and circles for η = 0.06.
Note the linear regime at finite η is delayed until larger F, giving way to the 3/2
regime prior to it.
less reliably. The reliability of the approach is hampered at larger η by the strong
divergence of trajectories explained in the previous section.
The pattern of escape trajectories for various trajectory parameters is shown in
Fig. 3.3 (see also appendix C for a survey of other interesting features of trajectories).
We plot several (∼ 10) optimal escape trajectories for various values of F and η as
thin dotted lines (see footnote 4 of Chapter 3 to recall the two meanings of the word
“escape”; here the first meaning is implied). Thick dotted lines denote the separatrix.
In each case, only one of the trajectories leads to the saddle.
A crossover between the ξ = 3/2 and ξ = 1 regime is defined to be such x at
which fits to the respective regions intersect. A crossover between the ξ = 1 and a
nonpower law regime was defined to be the point at which the fit to ξ = 1 regime
came closest to the ln (Q∗ /η) versus ln (x), see Fig. 3.1. Collecting data for S ∗ versus
x for various η we were able to map out the entire bistability region for locations of

Figure 3.3: (a),(b) Up → down escape trajectories for η = 0.001; (a) at 7% of the
hysteresis (i.e. x = 0.07), (b) at 0.001% of the hysteresis. (c),(d) Up → down escape
trajectories for η = 0.1; (c) at 80% of the hysteresis, (d) at 0.001% of the hysteresis.

Figure 3.4: Regions of different scaling behaviors for up → down transitions obtained
by scanning F at fixed η for multiple values of η. The meaning of the three scaling
regions is described in the text. The dashed line is a simple prediction of the crossover
between the ξ = 1 and ξ = 3/2 regime expressed by Eq. (3.3).
various scaling regimes and crossovers between them. The computations were made
down to η = 10−3 . The result of this work is depicted in Fig. 3.4. Three distinct
regions have been identified; region I: nonpower-law regime, region II: ξ = 1 regime,
and region III: ξ = 3/2 regime. We remind the reader that in Fig. 3.4 for up → down
transitions and in Fig. 3.9 for down → up transitions, the mapping was performed by
scanning F while holding η fixed. Note in Fig. 3.4 that as η → 0, the linear scaling
is found between x = 0 and some finite x. At small, but nonzero η, the region of
ξ = 3/2 subsequently grows. These results are in accord with the theorem of DSS
which states that close to a bifurcation which at zero η becomes nonlocal (i.e., the
low F bifurcation in our system), the underdamped regime (regime B) at nonzero η
is expected to give ξ = 1. Since the underdamped portion of the hysteresis extends
all the way to the bifurcation as η → 0, we expect the linear scaling to extend all the
way to the bifurcation in this limit.4

In practice we could go down to η = 0.001. In principle, strictly at η = 0 escape trajectories do
not exist; they only have meaning in the limit η → 0.

Figure 3.5: S ∗ /η versus x and η up to η = 0.175 which is about 30% of the hysteresis.
Fig. 3.1 are two slices of this surface.
The boundaries between the different scaling regions I, II and III are purely conventional - they have been defined in some convenient way, but do not correspond
to anything specifically physical. To emphasize this point, it helps to construct a
3-dimensional plot of S ∗ /η versus x and η - Fig. 3.5. This surface is smooth and
large parts of it happen to fit well to power laws. The regions when this is possible
appear to be much bigger then the regions at which the ξ = 1 and ξ = 3/2 scaling
were supposed to be applicable on the basis of approximations discussed in Section
2.2. A natural question is why this is so? To shed more light on this issue, let us
digress on the crossover from the ξ = 3/2 to ξ = 1 in more detail.
3.2.1.2

Discussion of the ξ = 3/2 to ξ = 1 Crossover

We now discuss the various ways in which the crossover from the ξ = 3/2 to the ξ = 1
regime may be defined. Notice that in Fig. 3.1 such crossover was defined as the point
of intersection of the fits to the respective regimes. Thus, to predict this crossover we

simply equate the asymptotic in the ξ = 3/2 regime, R = 4 3 2 ηδ 3/2 to the asymptotic

43
in the ξ = 1 regime R = 2ηδ (this is most accurate as η → 0). The result is
δ = η,

(3.2)

F = FBl + η,

(3.3)

which corresponds to

(FBl is itself a function of η). The crossover predicted this way is depicted by the
dashed line in Fig. 3.9. The fact that this method of predicting a crossover works well
tells us that the asymptotic expressions are correct. It does not teach us about the
crossover region itself: if it is a very wide region for example, the crossover defined
by equating the fits or predicted by equating the asymptotics is not a particularly
physically meaningful quantity. The width of the crossover region, for example, is a
more physically interesting quantity. Alternatively, we may also define the following
two crossovers: a “lower crossover” (LC) defined as a characteristic value of F (or δ
or x) above which the ξ = 3/2 regime begins to fail and an “upper crossover” (UC)
defined as a characteristic value of F (or δ or x) below which the ξ = 1 regime begins
to fail. The width of the crossover region is the difference between these two values.
Physically, we expect the LC to take place where the soft mode picture fails.
Therefore, a good characteristic estimate of this will be such δ at which the real
and imaginary parts of eigenvalues become equal. This can be shown to happen for
δ = η 3 /2, suggesting the cubic scaling. To predict the LC more accurately we turn
back to the type of analysis described in section 2.2.1. Let us rewrite Eqs. (2.56)(2.57) for the nonadiabatic part of the fast variable A = f − aB s. Near the low-field
bifurcation, these will have the form
Ȧ = −aB δ − 2ηA + a1 s2 + a2 As + a3 A2 + a4 s3 + a5 s2 A + a6 sA2 + a7 A3 + nx − aB ny ,
(3.4)
ṡ = δ + b1 s2 + b2 As + b3 A2 + b4 s3 + b5 As2 + b6 A2 s + b7 A3 + ny .

(3.5)

All coefficients ai and bi are functions of η. The position of the fixed points sets the
characteristic scaling of A to be δ and of s to be δ 1/2 (a generic feature of saddle-node

44
bifurcations), so we define new variables: s = δ 1/2 s0 and A = δA0 . Moreover, the
η-dependencies of the coefficients ai and bi near the low-field bifurcation are such
that when η is small compared to 1, further rescaling of s0 and A0 as s0 = η 1/2 σ
and A0 = η −2 φ yields a single dimensionless parameter µ = δ/η 3 ; such reduction is
approximate, but becomes asymptotically accurate for smaller η.5 In terms of µ, σ
and φ we have
dφ
dτ
dσ
dτ

φ2
+ µ3/2 3σφ2 + µ2 φ3 + Nφ ,
−1 − 2φ + σ 2 + µ1/2 σφ + σ 3 + µ 3φσ 2 −
(3.6)
3 2
= µ1/2 1 −
+ µ σφ − σ 3 + µ3/2
φ − 3σ 2 φ − µ2 3σφ2 − µ5/2 φ3 + Nσ ,
(3.7)

where τ = ηT . With the help of the definition of σ and Eq. (2.58) we note that the
zeroth-order part of the fast variable aB s is given by δ 1/2 η −1/2 σ and the correction due
to nonadiabaticity of the fast variable is A = δη −2 φ. We also see from Eqs. (3.6)-(3.7)
that when µ ¿ 1, both σ and φ are O(1). Hence comparison of the two terms says
that when µ ¿ 1, nonadiabatic part/adiabatic part ∼ (δη −2 )/(δ 1/2 η −1/2 ) = µ1/2 → 0
as µ → 0. In other words, in this regime the nonadiabatic part becomes unimportant
and the problem reduces to 1 dimension defined by the noisy dynamics of the slow
variable. To do this reduction explicitly, we first notice from Eqs. (3.6)-(3.7) that
at lowest order, the variable σ can be considered constant while φ(σ) ≈ 43 σ 2 − 12 +
nonadiabatic corrections. This φ(σ) is then substituted into Eq. (3.7). The resulting
expression can only be considered up to O(µ) since the higher-order terms must also
include the nonadiabatic corrections to φ to be complete [note that if the µ1/2 -terms
in Eq. (3.7) contained any φ terms, the complete expression within the adiabatic
approximation would only extend to O(µ1/2 ), not O(µ)]. The barrier height in the
From our definitions, ωF = ω0 1 + Qη
. Recall that in order for equations (2.36)-(2.37) to
apply, we must be working in the regime when ωF = ω0 (1 + small number). Therefore, if the limit
η → 0 is taken, it is implied that the limit Q → ∞ is taken first. Moreover, wherever the limit δ → 0
is taken, it is taken last, because we talk about approaching the bifurcation at constant η. However,
since η is always less then ηc ∼ O(1), the η is already small (compared to 1) and so wherever the
statements about “small” η are made, they do not imply the formal η → 0 limit.

45
resultant 1-dimensional potential is
4 2 1/2 √ 3/2
∆U =
µ + 2µ .

(3.8)

Again, the correction is within the adiabatic approximation; a nonadiabatic correction
would be higher order in µ. Taking into account that the Nσ has a diffusion constant
Dσ = δηD2 , we have

4 2 1/2 3/2 √ δ 5/2
R=
η δ + 2 5/2

(3.9)

[compare with Eq. (2.63)]. We see again that the correction term grows to be dominant when µ ∼ 1 (i.e., when δ ∼ η 3 ), but this is precisely when the nonadiabatic
part becomes comparable to the adiabatic part, and it turns out to be meaningless
to treat the problem as effectively 1-dimensional. Our analysis taught us that in the
regime when δ ¿ η 3 , the activation barrier R will scale as δ 3/2 to a good approximation. Outside the regime of δ ¿ η 3 the current analysis only allows us to conclude
that if ξ remains to be 3/2, the physics behind this scaling is not the physics of a
1-dimensional overdamped soft mode: the escape problem is not 1-dimensional there.
To answer the question of what happens outside of the δ À η 3 regime we chose to
compute the most likely escape path numerically and to calculate S ∗ on that path.

We looked at the difference between the zeroth-order theory, R0 = 4 3 2 η 1/2 δ 3/2 and
the numerically computed S ∗ . An example of this for η = 0.01 is shown in Fig. 3.6,
but to be consistent with out plots, we display R1 /η ≡ (S ∗ − R0 )/η versus x (i.e.,
S ∗ (δ) and R0 (δ) were reexpressed to be S ∗ (x) and R0 (x) and both were divided by
a constant η = 0.01, which of course does not affect the slope of the fit versus x).
This procedure was repeated for several values of η and in each case, a fit to the
function ∝ xp was made at low values of x. The first few points at lowest x where
the scatter was large were excluded; although the numerical calculations of S ∗ do not
have an inherent randomness, the data exhibits a scatter at low x due to the difficulty
of hitting the saddle and the discreteness of the time steps (the trajectory accelerates
exponentially in the initial stages). As before, we did not use a systematic criterion

Figure 3.6: R1 (x)/η ≡ (S ∗ (x) − R0 (x))/η plotted for η = 0.01. The straight line
indicates a fit within the range spanning from x ≈ 0.0013 to x ≈ 0.019; the slope of
the fit for this particular value of η is ≈ 2.51.

Figure 3.7: Main plot: the exponent p obtained from the fit of R1 (x)/η, plotted versus
η. Insert, the position of LC defined as that x at which 0.2R0 (x)/η = R1 (x)/η.

47
for the high end of the range of x over which the fit was made in accordance with the
assumption that the value of the fitted slope approaches an asymptotic figure when
the high end of the fitting range is not large. Values of p found this way are shown
in Fig. 3.7. As η → 0, p seems to approach ≈ 2.2, not far from p = 5/2 predicted
above. However, the crossover LC defined as that x at which 0.2R0 (x)/η = R1 (x)/η
scales linearly with η at small η, not cubically, as shown in the insert of Fig. 3.7.
The question of why the ξ = 3/2 scaling holds over a much larger region of
parameter space then expected based on the 1-dimensional soft mode picture remains
unanswered, and forms a good challenge problem for future work. One potential
hypothesis is that in the region where the 1-dimensional soft mode picture does not
hold true, the action along the MPEP actually grows only over a small portion of the
MPEP where the motion is essentially 1-dimensional. We provide a plot of the action
versus the distance along the trajectory, S(`) and the derivative dS/d` for a particular
value of parameters deep in region III, Fig. 3.8. The derivative plot clearly shows
that the action grows fastest in the last winding of the trajectory, but this is not
sufficient to conclude that the dynamics in that part of the trajectory can somehow
be described as effectively 1-dimensional.
Analogously to the lower crossover (LC), the upper crossover (UC) can be defined
as a characteristic value of the driving field below which the linear scaling begins
to break down. To predict this point, we must again compare the lowest and one
higher order terms in the theory for R(δ) or R(x) that applies in the ξ = 1 regime.
The effect of a nonzero η was clearly demonstrated in Fig. 3.2. The lowest-order in
dissipation (quasi-Hamiltonian) theory, Eq. (2.69) indeed predicts linear scaling for
low values of δ (or x), which according to Eq. (2.70) is given by 2ηδ. The region of
linear scaling persists at larger η, but due to the existence of the 3/2-scaling region,
this linear part has been pushed over to larger δ, and hence R takes lower values then
those predicted by the zeroth-order theory (see Fig. 3.2). Chinarov, Dykman, and
Smelyanskiy [81] worked out dissipative corrections to the quasi-Hamiltonian theory
to next order in η. Their result, translated to our notation reads: R = η (2δ − πη)
(compare with Eq. 2.70). The factor πη 2 is just that correction which takes this

Figure 3.8: (a) S(l) and its derivative for particular parameters in region III: here
η = 0.15 and F = 10% of the hysteresis, corresponding to µ ≈ 5; the soft mode
picture is expected to hold only for µ ¿ 1. (b) The plot of the MPEP along which
S(l) in (a) is computed; dashed curve denotes separatrix.
lowering into account. We expect the linear scaling to start breaking down when the
η for
correction reaches some fraction r of the zeroth-order term, i.e., when δU C = 2r

low η.
To extract δU C from the data we followed the same procedure used to extract
δLC : we subtracted the zeroth-order part R0 (x) from S ∗ (x)/η. However, recall what
S ∗ (x) does: it has an initial regime with the leading exponent 3/2, followed by a
linear regime for higher x, followed by a nonpower law regime for higher x still. On
the other hand, R0 starts out linear: it does not have the 3/2 for the up → down
transitions (for down → up transitions, even R0 has a 3/2 regime). Therefore, the
function (S ∗ (x) − R0 (x))/η grows for low x, slows down for intermediate x, and
increases for yet larger x past the linear regime. Note that in the intermediate regime

49
the function S ∗ (x)/η continues to grow because the exact slope in the linear regime
is slightly different from the zeroth-order result! Thus, because (S ∗ (x) − R0 (x))/η
does not plateau in the intermediate, linear regime, it is hard to make quantitative
predictions. We attempted to extract an average value in the intermediate regime,
although is is somewhat ambiguous because it involves a conventional choice of the
interval over which this averaging is performed. Nevertheless, equating this correction
to the zeroth-order term yields δU C which does depend linearly on η.

3.2.2

Down → Up Transitions

The scaling regimes for the down → up transitions were also mapped out over the
entire bistability region. The product of this mapping is shown in Fig. 3.9. The
computations were made down to η = 0.005.

Figure 3.9: Regions of different scaling behaviors for down → up transitions obtained
by scanning F at fixed η for multiple values of η. Region I corresponds to a nonpowerlaw regime. Region II corresponds to the ξ ≈ 1.3 regime. Region III corresponds to
the ξ = 3/2 regime. Fuzzy regions correspond to parts of the (F, η) plane which
proved very difficult to map reliably due to the high divergence of trajectories.
A portion of the surface S ∗ /η(x, η) is presented in Fig. 3.11 and two cuts at a fixed η

Figure 3.10: Main plots: S ∗ (x)/η (dots) and fits (lines) for down → up transitions.
Fits to the ξ = 3/2 and ξ ≈ 1.3 regimes was made by analyzing ln (S ∗ /η) versus
ln x. In the region where the simple shooting method was unreliable for hitting the
saddle due to the high divergence of trajectories, a bracketing method was used. The
results from using this method are denoted by squares in plot (b). Inserts: S ∗ (x)/η
(dots) and theory (lines) based on the quasi-Hamiltonian theory, Eq. (2.69). Note:
this theory predicts the ξ = 3/2 scaling, as was first pointed out by [77]; in the light
of the recent work by [21] this was related to the nonlocality of the bifurcation at zero
η.

Figure 3.11: S ∗ /η versus x and η up to η = 0.22 which is about 33% of the hysteresis.
Fig. 3.1 are two slices of this surface.
appear in Fig. 3.10. The low amplitude basin does not exhibit the nonlocal geometry
discussed in Sec. 2.3 upon approach to the bifurcation at zero η, and correspondingly,
the linear scaling is not found in the down → up transitions at any η (see Fig. 3.10),
which is in accord with the DSS theorem. Instead, the ξ = 3/2 scaling gives way to a
different scaling relationship further from the bifurcation, with an exponent ξ ≈ 1.3.
We currently do not have a physical explanation for the origin of such exponent; there
may not be such an explanation, and 1.3 may simply be a power-law interpolation of
1 or 2 decades of numerical data. Note that region III in Fig. 3.9 (the ξ = 3/2 region)
extends all the way to zero η - this is unlike the case of up → down transitions, where
the ξ = 3/2 region disappears as η → 0. These results are also in accord with the
theorem of DSS which states that close to a bifurcation which is local at η = 0, the ξ is
expected to be 3/2 in the underdamped regime at nonzero η. Since the underdamped
portion of the hysteresis extends all the way to the bifurcation as η → 0, we expect the
3/2 scaling to extend all the way to the bifurcation in this limit. At finite η, the region
of the overdamped 3/2 scaling and this underdamped 3/2 scaling merge smoothly into
each other. Had we assumed that ξ = 3/2 is due only to the overdamped soft mode,

52
we would follow the nearly adiabatic analysis of the previous section, and we would
find that for low η there exists another dimensionless parameter ν = δ 1/2 η −2 . When
ν ¿ 1 such analysis predicts ξ = 3/2. When ν ∼ 1, the reduction to 1 dimension due
to the adiabaticity becomes impossible and predictions based on this method can not
be made. The numerical observation that ξ remains 3/2 beyond ν ¿ 1 confirms that
the 3/2 scaling is not a consequence of reduction to an overdamped 1-dimensional
soft mode.
In the previous section we offered three definitions of the 3/2-to-1 crossover. In
the case of up → down transitions, we do not have a theory for the ξ ≈ 1.3 scaling
regime, so the crossover from the 3/2 scaling can be defined here only as a point at
which the 3/2 scaling begins to break. Since the 3/2 scaling is predicted by the quasiHamiltonian theory, we look at next order corrections due to Chinarov, Dykman,
and Smelyanskiy [81]. Translating the results of that work to our notation says that
ηδ 3/2 + 1.06 31/4
η 3 |ln η| δ 3/2 (compare with Eq. 2.70).
close to the bifurcation R ≈ 31/4

Evidently, the correction is to simply shift R(δ) by a constant offset (see Fig. 3.10
(b) where this effect is clearly pronounced throughout the entire R(δ)). Because the
correction is so small for small η, its effect does not become significant until a rather
large η. Equating the zeroth-order term to the correction term predicts this value
of η to be approximately 1.5 > ηc . Numerically, this delayed effect is observed (see
Fig. 3.9): the crossover does not seem to depend on η significantly until η ≈ 0.4ηc .
However, we state this with caution because for η larger than ≈ 0.4ηc we were not
able to extract the position of the crossover accurately due to the high divergence of
trajectories (the fuzzy region in Fig. 3.9).

Chapter 4
Consequence of Breaking of
Detailed Balance
In this chapter we focus on the probability density function (PDF) in the (Q, P )
phase space. We first bring to light the notion that a steady-state nonequilibrium
system that satisfies detailed balance behaves essentially like an equilibrium system.
It is the breaking of detailed balance that causes drastic differences to take place. We
show that in our oscillator system, the detailed balance is broken at any finite value
of the driving F, so in this system the lack of detailed balance is equivalent to the
lack of equilibrium; in other words the detailed balance is fragile in our system. We
then present our results on the formation of nondifferentiable features in the PDF in
the limit of vanishing noise strength1 and find that the smoothness of the the PDF is
also fragile. Section 4.1 introduces the meaning of detailed balance for Fokker-Planck
models, states the consequences of detailed balance, proves an interesting relationship
between the fluctuational and relaxational trajectories in the presence of detailed
balance and finally lists the possibilities that might occur if detailed balance does not
hold. Section 4.1 may seem somewhat formal, but the necessary information related
to subsequent sections is given in its last several paragraphs. The section that follows,
Sec. 4.2, studies the listed consequence of breaking of detailed balance specifically for
our nonlinear oscillator.
The resulting physics of is qualitatively analogous to the zero wavelength limit in optics (ray
optics) or to the quantum tunneling in limit of ~ → 0.

4.1

Detailed Balance

There is a long history of research aimed at generalizing the methods of equilibrium
statistical physics. One of the central questions is the feasibility of a nonequilibrium
system to attain a Boltzmann-like distribution2
ρ(r) = Ae−φ(r) .

(4.1)

Depending on the problem, φ may also depend on some measure of effective temperature or noise strength. Narrowing the discussion to Fokker-Planck models, the
necessary and sufficient condition for the existence of such a φ appears to be the
presence of the detailed balance with respect to a time reversal transformation [29],
as was shown by Graham and Haken (GH) in 1971.
Consider a 2-point distribution function ρ(r2 , t2 , r1 , t1 ). Under a time-reversal
transformation, the variables r transform in some way:
t → t̃ = −t

(4.2)

r → r̃

(4.3)

For example, when r = (x, p), i.e. position and momentum, then x̃ = x and p̃ = −p.
Detailed balance means the following [2]:
ρ(r2 , t2 , r1 , t1 ) = ρ(r̃1 , t2 , r̃2 , t1 ).

(4.4)

For physical processes such as atomic collisions for example, the criteria for this
condition come from the fundamental equations of motion. For us, the discussion is
restricted to the cases when ρ evolves according to FPE (the two-point distribution
function can be written as a product of a one-point distribution and a conditional
transition probability distribution, both of which obey FPE). Consider a FPE of the
form Eq. (2.41) with D replaced by a spatially dependent diffusion matrix D (assumed

Due to the abundance of vectors and indices, vectors will be represented by boldface characters
in this chapter.

55
symmetric) with components Dij (GH also factor out 1/2 from it). According to GH,
if detailed balance holds and we require the solution to the stationary FPE of the
form e−φ(r) , then certain potential conditions will hold. To state these conditions, it
helps to introduce a few definitions. First, let
K̃i (r) = ²i Ki (r̃(r)),

(4.5)

D̃ij (r) = ²i ²j Dij (r̃(r)),
 +1 if r̃ = r ,
²i =
 −1 if r̃ = −r

(4.6)

(4.7)

(no summation over indices is implied). The K̃ and D̃ represent transformations to
K and D under time reversal, just as r̃ represents transformation to r under time
reversal. Moreover, we may break the K into the part that changes sign under the
transformation (the reversible part, R) and the part that does not change sign (the
irreversible part, I): K = R + I. Then
(Ki (r) + K̃i (r)),
(Ki (r) − K̃i (r)),

Ii =

(4.8)

(4.9)

Given these definitions, the first potential condition is a condition of the diffusion
matrix and reads
Dik (r) = D̃ik (r).

(4.10)

The second condition dictated by the detailed balance says that the irreversible part
of the probability current must be zero:
Ii ρ(r) −

(Dik ρ(r)) = 0
∂rk

(4.11)

When ρ = Ae−φ is substituted into this equation, we get the second potential condition:
∂φ
= Dli−1
∂rl

∂Dik
− Ii
∂rk

(4.12)

56
This serves as an equation for a smooth φ.3 Finally, the third potential condition
reads
0=

∂Ri
∂φ
− Ri
∂ri
∂ri

(4.13)

With the help or Eq. (4.1), this can be recast as
∂(Ri ρ)
= 0,
∂ri

(4.14)

which says that the reversible part of the drift must support the time-independent
ρ(r). All together these are conditions on the diffusion matrix and the drift field.
The converse also holds true: if the potential conditions hold, then there is detailed
balance and there exists a smooth φ obtained from Eq. (4.12) by quadrature. This
is remarkable given that the system may actually be driven away from equilibrium!4
We now derive an interesting consequence of detailed balance. To do this, it is
helpful to first note the effect of generalizing a scalar D to a matrix on the Classical
Mechanics described in Sections 2.1.3 and 2.1.2. Let D = Dd - we factored out the
noise strength D, just a scalar control parameter. Then the Lagrangian becomes
L = d−1
µν (ṙ − K (r)) (ṙ − K (r))

= (ṙ − K)T d−1 (ṙ − K),

(4.15)
(4.16)

(now sum over repeated indices), the Hamiltonian in Eq. (2.48) generalizes to
H(r, p) = dµν pµ pν + Kµ (r)pµ
= pT dp + K(r)p

(4.17)
(4.18)

(recall that r stands for the vector (Q, P ) in our problem). In the case of a unit
matrix d we recover the old expressions for L and H. The auxiliary equations of
From Eq. (4.12) it follows that ∂r∂φ
= ∂r∂φ
i ∂rj
j ∂ri
The potential conditions were known before GH, but it was their work [29] that explained the
physical meaning behind these conditions to be detailed balance.

57
motion read
∂H
= 2dkµ pµ + Kk (r),
∂p
∂H
∂dµν µ ν ∂K µ
= − k =−
p p −
pµ ,
∂r
∂rk
∂qk

ṙk =

(4.19)

ṗk

(4.20)

and finally the momenta are related to the action S in the usual way: p = ∇S. The
distribution is still given by
S(r)
ρ(r) ∼ exp −

(4.21)

We compare this to the definition of φ in Eq. (4.1) and conclude that S = Dφ, so ∇S =
D∇φ. We now specialize to the case of spatially independent d. Then Eq. (4.12) tells
us that p(r) = −d−1 I(r). We substitute this into EOM for ṙ, Eq. (4.19), and obtain
ṙk = −2Ik (r) + Kk (r).

(4.22)

ṙF L = −K̃(r).

(4.23)

In other words,

The subscript “FL” is to remind that this equation describes a fluctuational trajectory. On the other hand, the relaxational trajectory obeys
ṙREL = K(r).

(4.24)

Therefore, the flow field that describes the fluctuational trajectory (or escape trajectory) is the negative of time-reversed flow field that describes the relaxational trajectory! In our particular system, Eqs. (2.36 - 2.40), the diffusion matrix is Dij = Dδij
with constant D, and in the next section we will see explicitly how this conclusion
holds true by calculating analytically the relaxational and the fluctuational trajectories when detailed balance holds.

The next logical question is, what are the consequences of violation of detailed
balance? Specifically, what happens to the PDF ρ(r)? It is reasonable that in general
ρ(r) can still be written down in a Boltzmann-like form; we see this from the pathintegral solution to the FPE: unless the fluctuations around the optimal (classical)
path contribute to the exponent, the distribution will still have the form ρ(r) =
ρ0 (r) exp [−S(r)/D]. However, the prefactor may now depend on r and it may even
be singular [82] and the action S may in general be a multivalued function of r if
there exists more then one escape trajectory reaching this r from r0 . Generically, these
trajectories will yield a different value of the action S1 and S2 , and in the limit of small
noise exp [−S1 /D] may drastically differ from exp [−S2 /D] although the potential
becomes smoothed out at a finite noise strength [30]. Thus, in the limit of vanishing
noise strength one must always choose that trajectory which leads to the least value
of the action, a procedure that may lead to a nonsmoothness (or nondifferentiability)
of φ with respect to r. This effect was described by several authors; Graham and Tél
in [30] appear to be one of the earlier authors to describe this particular mechanism
of formation of nondifferentiability of φ. They point out a crucial geometric fact of
the auxiliary dynamical system: in order for a projection of two trajectories in the
r − p phase space to meet at one point in the r space, the Lagrangian manifolds on
which these trajectories live must “fold over”: the trajectories cannot intersect on the
manifold, in virtue of the existence and uniqueness of solutions of dynamical systems,
so the only way for their projections unto the r plane to intersect is for the manifold to
fold over. Dykman, Millonas and Smelyanskiy further clarify the geometric details of
this in [31] (one is highly encouraged to have a look at Fig. 2 in this work for a cartoon
explanation of the relationship between folds, caustics, switching lines and trajectories
reflecting off of caustics). Relating to Catastrophe theoretic classification of folds in
Lagrangian Manifolds (Appendix 12 of [76] and [83, 84]), they note that such foliation
generically happens in a cusp manner, also known as Whitney’s tuck. Because the
intersections of neighboring trajectories happen at the folds, the projections of the
folds onto the r space map out a pair of caustics (see [31] or Maier and Stein’s

59
paper [82]).5 Dykman, Millonas and Smelyanskiy also describe a generic form of the
caustic near a saddle fixed point of the deterministic system (i.e. of K). Although the
caustics are not experimentally observable, the “switching line,” which is a curve at
which S1 = S2 [30, 31, 85, 86] is the feature which should be observable in experiments
(see references in [31]). The nondifferentiability of S(r) takes place along switching
lines.
In the 1980s Graham and Tél undertook the study of caustics in various models
and made an attempt to place conditions on the structure of the FPE which would
entail smooth potentials (assuming detailed balance is broken). Their search for such
conditions centered on the search for integrability criteria [32], although they state
that integrability is not a sufficient condition, because it does not prevent multivaluedness. Yet, integrability in their point of view was sufficient to prevent infinite
foliation of Lagrangian manifolds, a phenomenon also known as wild separatricies
[30, 87]. A general condition that ensures integrability of the auxiliary system associated with the FPE Hamiltonian H appears elusive. In the words of Graham and Tél
[30], “methods testing for integrability on certain hypersurfaces in phase space, like
E = 0, would be better suited for the problem at hand, then tests for integrability in
the entire phase space, but methods of this type remain to be developed”.
Most of the efforts to describe caustics appear to rely upon the existence of a
fixed point or a limit cycle. Among more recent works, the formation of caustics in
the presence of a saddle fixed point has been investigated by Maier and Stein in [82].
This work studies conditions when the caustic is attached to the saddle and when it
breaks off as a control parameter in the model is varied. The case when the system
contains a limit cycle has also been understood in great detail by [88] for example.
All of Graham’s works cited here appear to focus on models with limit cycles as well.
To my knowledge, a work by Dykman, Millonas and Smelianskiy [31] appears to be
the only paper in which the situation without any limit sets (aside from the attractor
from which the trajectories escape) is considered. The system they describe happens
to be the Duffing oscillator in the monostable regime (they also consider the bistable

In mechanics or ray optics, a caustic is an envelope of intersecting trajectories or rays.

60
regime where there is a saddle; they study conditions for detachment of the caustic
from the saddle). Indeed, as described in Chapter 1, for F < Fl there exists only
one single attracting fixed point. They write down the normal form for ρ(r) in the
vicinity of the cusp point - the point where the pairs of caustics meet. What is not
discussed in this work is the evolution of the caustic system as the driving is turned
up, and most importantly, the mechanism of the formation of the caustic pair. The
questions are, is there a finite threshold F at which the caustic appear and what
serves to nucleate the caustics as F is turned up. We answer these questions here.
We make several surprising discoveries: (i) the caustics do exist; they appear
in pairs as expected on the grounds of general catastrophe theory [31, 83, 84] and
appendix 12 of [76]; (ii) there does not appear to be a threshold F for the appearance
of a caustic; (iii) the radial distance R (not to be confused with R above) to the cusp
point as a function of F appears to be a power law: R ∼ F −p . Thus, as F is turned
on from zero the caustic appears at infinity and is brought to smaller radii for larger
F. We now present these findings, which are mostly numerical, and describe our (so
far largely unfruitful) attempts to explain them analytically.

4.2

Fragility of Detailed Balance and Formation
of Caustics in the Monostable Regime of the
Oscillator

Before moving on, it is instructive to check if and when our system represented by
Eqs. (2.36 -2.40) satisfies detailed balance [i.e., satisfies the potential conditions and
admits a smooth φ(r) or S(r)]. First, note that the variables of the Duffing oscillator
have the following obvious time-reversal properties: x̃ = x and x̃˙ = −x. According
to the definition of Q and P , Eqs. (2.28)-(2.29), these will have the following timereversal properties: P̃ = −P and Q̃ = Q. Then, using the definition of KQ and KP ,
Eqs. (2.26)-(2.27) and the definitions in Eqs. (4.5)-(4.8) we get

IQ = −ηQ,

(4.25)

IP = −ηP,

(4.26)

RQ = −P + P (Q2 + P 2 ),

(4.27)

RP = Q − Q(Q2 + P 2 ) + F.

(4.28)

Given this, we construct the potential conditions. The first potential condition,
Eq. (4.10), is satisfied automatically. The second potential condition, Eq. (4.12)
says
∂φ
ηQ
= − IQ =
∂Q
∂φ
ηP
= − IP =
∂P

(4.29)
(4.30)

This is indeed an equation for φ with the simple solution
φ=

(Q2 + P 2 ).
2D

(4.31)

And the third potential condition expressed in Eq. (4.13) says
∂RQ
∂φ
ηP Q £
− RQ
= 2P Q +
1 − (Q2 + P 2 ) ,
∂Q
∂Q
¤ ηFP
∂RP
∂φ
ηP Q £
0 =
− RP
= −2P Q −
1 − (Q2 + P 2 ) −
∂P
∂P
0 =

(4.32)
(4.33)

Adding the two equations we get
ηFP
= 0.

(4.34)

Thus, detailed balance holds only when F = 0, in other words the detailed balance
is fragile!
The auxiliary EOM were explicitly stated in Sec. 2.1.3 [Eqs. (2.51)-(2.54)]. Some-

62
times, it will prove useful to work in polar coordinates, so in terms of the radial
coordinate r and the angular coordinate θ, the Lagrangian is
i F2
´2 ¸ F h³
rθ̇ − r + r3 cos θ + (ṙ + ηr) sin θ +
L=
(ṙ + ηr) + rθ̇ − r + r
the Hamiltonian is
H = p2r +

´ F2
p2θ
pθ
ηrp
(1
sin
cos
r2

(4.35)

and EOM are
ṙ = 2pr − ηr + F sin θ,
2pθ
θ̇ =
+ (1 − r2 ) + cos θ,
2p2θ
pθ
ṗr =
+ ηpr + 2rpθ + F 2 cos θ,
r ³
r´
pθ
ṗθ = −F pr cos θ − sin θ .

4.2.1

(4.36)
(4.37)
(4.38)
(4.39)

Exact Solution at Zero F

We can exploit the angular symmetry of the situation: when F = 0, the Hamiltonian
g is rotationally symmetric, g = 41 (Q2 + P 2 − 1)2 , so the 2D flow is symmetric and the
set of escape trajectories must also be symmetric!6 This is where polar coordinates
are most useful. Because of this symmetry, pθ = const: this is seen directly from
Eq. (4.39) and represents conservation of angular momentum, but this const must be
0 due to the fact that the trajectory emanates from the fixed point. Because of this,
and the fact that the motion takes place on the H = 0 surface, Eq. (4.35) reveals
∂S
that pr = 0 (relaxational paths) or pr = ηr (fluctuational paths). Now, ∂q
= pi , so
To remind the reader, this refers to the most likely trajectories along which the fluctuation from
the attracting fixed point to a point (Q, P ) takes place.

the action is S(r) = η r2 = η Q +P
, and pQ = ηQ, pP = ηP . Therefore,
∂g
∂P
∂g
Ṗ = ηP −
∂Q

(4.40)

Q̇ = ηQ +

(4.41)

This is consistent the prediction we made earlier from general considerations in
Sec. 4.1. When our nonlinear oscillator satisfies detailed balance, the escape trajectories are obtained from relaxational trajectories with a reversed sign of friction
[indeed, this is K̃(r), Eq. (4.23), for this particular system]. We also learn that the
distribution is given (aside from a normalization factor) by
³ η
ρ(Q, P ) ∝ exp −
(Q2 + P 2 ) ,
2D

(4.42)

which we already derived (Eq. (4.31)) via the more general method.
For the record, I write out the zeroth-order (F = 0) solution to Eqs. (4.36)-(4.37)
explicitly because we will need it below. It satisfies ṙ = ηr and θ̇ = 1 − r2 (remember
that pθ = 0 for F = 0).
· µ ¶
r2 − r02
ln
θ(r) = θ0 +
r0

(4.43)

r(t) = r0 eηt ,

(4.44)

or in terms of time,7

θ(t) = θ0 + t −

r02 ¡ 2ηt
2η

−1 .

(4.45)

I switch notation from T to t for now, and it will be understood that t is not the original time
in the Duffing equation, but the slow time in the amplitude equations, formerly known as T .

64
Or
Q(t) = eηt [Q0 cos ψ(t) − P0 sin ψ(t)] ,

(4.46)

P (t) = eηt [Q0 sin ψ(t) + P0 cos ψ(t)] ,
Q2 + P02 £ 2ηt
e −1 .
where ψ(t) = t − 0
2η

(4.47)
(4.48)

The solution is designed to arrive to point (Q0 , P0 ) at time t = 0. We will want to
think of the potential φ as a function of the field point (Q0 , P0 ) and then drop the
subscript.

4.2.2

An Aside

We can actually extract the result given by Eq. 4.42 from more physical arguments
directly for the Duffing oscillator - this will underline the usefullness of potential
conditions discussed in prior sections. In the absence of driving, the Duffing oscillator
is descried by q̈ + 2Γq̇ + ω02 q + γq 3 = f (t), Eq. (2.1). Clearly, in thermodynamic
equilibrium at temperature T (not to be confused with the slow time), q and q̇ obey
the Boltzmann distribution
ρ(q, q̇) ∝ exp −

q̇ 2
+ ω02 q2 + γ q4

(4.49)

Now we can switch to P and Q via the transformation in Eqs. (2.28)-(2.29). The
q̇ 2
+ ω02 q2

part gives
8 ω03
(ωF − ω0 ) P 2 + Q2 + O
2 3η γ

2η

(4.50)

³ ´
where the O

2η

term comes form the fact that q̇ contains ωF in its definition

(Eq. (2.28)) while q 2 has an ω02 in front of it. The γq 4 term is of the form
η 2 ω02
(ωF − ω0 )2 × O(P 4 , Q4 , ...).
2 η γ

(4.51)

65
Thus (quartic term)/(quadratic term) is of order ωFγω−ω
= 2η² . Thus,

ρ(Q, P ) ∝ exp − 0 (Q2 + P 2 ) + O
2D

2η

¶¶

(4.52)

where
D0 =

3γ
× (2ΓkT ).
16ω0 (ωF − ω02 )

(4.53)

This equals D if B = 2ΓkT which is just an instance of a Fluctuation-Dissipation
³ ´
relationship. The O 2η² correction comes from the fact that the displacement is
assumed small and therefore q samples mostly the linear part of the potential. It
³ ´
may appear that because of this O 2η² , the ρ obtained with this method does not
³ ´
correspond exactly to the ρ obtained above, which did not mention O 2η² . However,
that derivation of ρ was based upon the slow amplitude equations, Eqs. (2.26)-(2.27)
which were obtained under the assumption that displacement q is small and samples
mostly the linear parts of the Duffing potential. In other words, in the amplitude
equations the ²-corrections have also been neglected.
What we have just seen is that if a Boltzmann distribution holds in the phase
space of the oscillator, then it also holds in the canonically transformed phase space.
From this point of view, it is clear why the new distribution will not remain given
£ η
(Q2 + P 2 ) : this is because the Boltzmann distribution in phase
by strictly exp − 2D
space does not strictly apply as soon as the time-dependent driving is introduced. By
definition, the problem becomes nonequilibrium since the notion of full exploration
of phase space does not apply: the system does not have enough time to explore the
phase space as the applied field changes it. But what is the picture from the point of
view of time-independent amplitude equations, Eqs. (2.26)-(2.27)? This is where the
ideas from Sec. 4.1 come in. Suppose we are just given the Fokker-Planck equation
without realizing that it came from a Duffing oscillator. The FPE is autonomous,
so we can not use time dependence to say that the Boltzmann distribution does not
apply, but we can use potential conditions as a “litmus test” for detailed balance and
hence for the possibility of a smooth φ(Q, P ). As we just saw, the potential conditions

66
fail at nonzero F. We will now see how φ(Q, P ) fails to remain smooth in practice.
The subject of Sec. 4.2.3 and Sec. 4.2.4 will appear in the upcoming publication [89].

4.2.3

Phenomenology of Caustics in the Monostable Regime

In this subsection we will study the pattern of escape trajectories at ever increasing
driving strength F. We recall that the exponential part of the distribution ρ(Q, P )
is given by the action of trajectories along which the system is most likely to escape
the attracting fixed point (or attractor) to the point (Q, P ). These are the “escape
trajectories” which form the subject of discussion in this section. The word “escape”
refers to the escape from the attractor to (Q, P ), not to be confused with the escape
from a basin which was relevant in Chapter 3. We also recall that one may refer to
these as “fluctuational” trajectories, while the trajectories along which the system
is most likely to return from (Q, P ), which are just the trajectories of the noise-free
dynamics of the oscillator, may be called “relaxational” trajectories.
Unless otherwise stated, the phenomenology will be demonstrated at fixed damping η = 0.4. When F = 0, the pattern is rotationally symmetric, Fig. 4.1. This
fact was proven above in terms of the symmetry of ρ(Q, P ), Eq. (4.42), and in terms
of the explicit analytical form of the trajectories, Eq. (4.43). All calculations of
trajectories in this section are terminated either when a stopping radius is reached
or when a trajectory reaches a caustic, as explained shortly. The zeroth-order solution, Eqs. (4.44)-(4.45) suggests that trajectories are accelerating exponentially, which
means that to provide an accurate solution valid at larger radii, the time steps must
decrease dramatically. The choice of the stopping radius was dictated mostly by the
position of interesting features (caustics). As we will see, these interesting features
recede to larger radii as F decreases, and so calculations at very low F required very
small time steps (this would be an ideal place for a variable time-stepping algorithm,
but we used a 4th order Runge-Kutta). For a given stopping radius, the time steps
were decreased until the results did not change significantly. Then more extensive
calculations were performed at this saturated value of time steps.

−1

−2

−3

−4
−4

−3

−2

−1

Figure 4.1: The pattern of 10 escape trajectories at F = 0. Each trajectory undergoes
approximately 3 revolutions after reversing its direction of rotation at the slow-down
radius of 1.
Fig. 4.2 displays the trajectory pattern at F = 0.005, clearly showing the breaking
of the rotational symmetry: the density of trajectories becomes nonuniform, and some
of them start to “jam” into each other. As F is increased further, the radius at which
this jamming takes place decreases. One way to find out where this happens exactly is
to consider pairs of trajectories with a similar neighboring parameter and to observe
the position of their intersection. There is a more elegant and powerful way. In our
discussion of why projections of trajectories intersect, we mentioned that the reason
for this phenomenon is the folding of Lagrangian manifolds in the (Q, P, pQ , pP ) space.
As two trajectories with an infinitesimally different trajectory parameter pass over
this fold, their projections will intersect at the fold; excellent cartoons of this can be
found in [31] and also in the beginning of [19]. The Lagrangian manifold is the surface

−1

−2

−3

−4
−4

−3

−2

−1

Figure 4.2: The pattern of 10 escape trajectories at F = 0.005.
p(r), i.e., (pQ (Q, P ), pP (Q, P )). The fold is characterized by such location where the
tangent plane to the surface is perpendicular to the (Q, P ) plane. This happens when
∂(Q, P )
= 0,
∂(pQ , pP )

(4.54)

in other words, when the 1-to-1 correspondence between ∆p and ∆r is lost. Alternatively, we can say that
J=

∂(pQ , pP )
= ∞.
∂(Q, P )

(4.55)

Thus, knowing p(r), the caustics can be mapped out at locations where J becomes
infinity. In general, the web of caustics may be very complicated, even of infinite
complexity in the case of “wild separatricies” [87, 88].8 The function p(r) is not

This type of behavior has also been found in our system in the bistable regime. It can be seen
in the pattern of trajectories that hit the separatrix: as narrower and narrower ranges of trajectory

69
known a priori, it is something that we are in fact seeking; only p along a given
trajectory can be calculated. However, as the equations of motion for pQ , pP , Q and
P are integrated in time, it is also possible to simultaneously integrate the equations

∂p

∂ S
of motion for second derivatives such as ∂Q∂P
= ∂PQ , etc. These are called the time-

dependent Riccati equations, which are derived and stated by Maier and Stein [74].
Ṡij = −

∂ 2H
∂ 2H
∂ 2H
∂2H
Sik Sjl −
Sik −
Sjk −
∂pk ∂pl
∂rj ∂pk
∂ri ∂pk
∂ri ∂rj

(4.56)

These equations are not only of practical use; they have a rather deep differentialgeometric interpretation [74], which is outside the scope of this thesis. With the help
of these and the other four EOM we can track the evolution of J: all together, there
are 7 coupled ODEs since the ODE for mixed second partial derivatives of S are the
same. The initial conditions are obtained from the knowledge of the surface p(r) in
the vicinity of the attractor, for example, from the knowledge of the eigenvectors.
We now proceed with the demonstration of the effects of increasing the strength
of driving. At F = 0.008 the region of increased density of trajectories move in a
little bit closer to the origin. This is shown in Fig. 4.3 along with the portrait of the
caustic. We see the cusp structure, which is expected on the generic basis (see [31],
the Appendix 12 of [76] and [83, 84]).
Another word of technicality is in order. Although caustics are projections of folds
in the Lagrangian manifolds unto the (Q, P ) plane, the dynamics on the manifolds
is obviously not trivial. Two trajectories with different shooting parameters may hit
a caustic after going through very different dynamics; it is not obvious what each
of these two trajectories do before hitting a caustic, if it is destined to do so. In
particular, a stopping radius which is not large enough may prevent some trajectories
from reaching a caustic. In Fig. 4.3 for example, the portion of the caustic shown
lies around the radius of 4, but in mapping out this caustic, we evolved trajectories
choosing the stopping radius to be 6.
As F is increased further, the cusp structure continues to pull in closer to the
parameters are considered, the pattern repeats itself. This fractal structure is believed to be caused
by the proliferation of folds, one on top of the other in the vicinity of the saddle.

−1

−2

−3

−4
−4

−3

−2

−1

Figure 4.3: The pattern of 10 escape trajectories at F = 0.008 and the pair of caustics
originating at the cusp.
origin. Some caustic and trajectory patterns at larger F will be shown below. For
now we summarize the relationship between the radius at which the cusp is located
versus F for two values of η. This is shown in Fig. 4.4.
Two power laws R ∼ F p with a crossover appear. The nature of the crossover
is not obvious, but the discussion below may hint that it may be related to the
appearance of a different caustic around the attractor. For η = 0.4, p ≈ −0.502 at
very low F and p ≈ −0.386 at larger F. The second power law appears to hold rather
well for F as large as 0.39 (not shown).9 For η = 0.1, the two values of p ≈ −0.456
at very low F and p ≈ −0.242 at larger F. We also observe that the power laws at
lower F are less sensitive to η then the power laws that hold after the crossover, at

If the distance is measured not from the origin, but from the attractor, which does not sit
precisely at the origin, the points at larger F fall off the linear trend.

Radius

-6

-4

-2

Figure 4.4: Log-log plot representing the radius of the cusp, measured from the origin
versus the strength of the driving F for η = 0.4 circles and η = 0.1 diamonds. The
solid and dashed lines represent fits to the two power-law regimes.
larger F.
As F increases, the caustics also evolve to become more complex. Fig. 4.5 conveys
a sense of this evolution with increasing F. In each case, the caustics were mapped
out by evolving 5000 trajectories with the stoping radius of 6. At F = 0.05 we notice
the appearance of the second caustic around the attractor. Its tight spiral structure
is seen in greater more detail in panel (b). There are also occasional loose red dots,
indicating that in addition to the familiar cusp-type caustic and the new spiral caustic
around the origin, there are other, sparsely visited pieces of caustics. At higher F the
spiral caustic around the attractor appears to be the continuation of the inner branch
of the familiar pair of caustics, panels (c) and (d). Finally, the last two panels, (e)
and (f), depict the situation for F = 0.395, which correspond to the bistable regime:

1.5

(a)

(c)

0.5

−1

−0.5

−2

−1

−3

−1.5

−4

−2

−1
−2

(b)

(e)

−1

−2

(d)

(f)

0.2

0.1

0.5

−0.1

−1

−0.2
−0.2

0.2

−0.5
−0.5

0.5

−2

−1

Figure 4.5: Caustics systems at various F. Panels (a) and (b): F = 0.05, panels (c)
and (d): F = 0.2, and panels (e) and (f): F = .395, all at η = 0.4. In panel (f) we
also show a set of escape trajectories from one of the attracting fixed points, one for
a different trajectory parameter. The fixed points and the separatrix are also shown,
with the saddle laying on the separatrix.
two more fixed points, a saddle (in the middle) and another attractor appear. We
observe that the familiar caustic does not seem to be connected to the saddle.10 The
phenomenology of the caustic structure and of the pattern of trajectories can be come
quite complex at larger F. One may ask, for example, what types of trajectories lead
to the circular caustic around the origin and whether these ever interact with the
cusp-type caustic that we discussed. Indeed, there is plenty of room for future work
to study this phenomenology. On the other hand, our goal here was to understand
the features which are expected to be generic, such as the scaling properties of the
caustics as F → 0 and to answer the question whether there is a finite threshold F
10

Indeed, the ratio of repulsive egenvalues of the saddle is ≈ 9.9 here, so in light of [31], the caustic
is not expected to be attached to the saddle for this particular choice of η and F.

73
at which the caustics appear. This is the main thrust of numerical Sec. 4.2.3 and the
analytical Sec. 4.2.4 that follows.
Varying η at fixed F also reveals a scaling which appears to develop asymptotically
in η. This is shown in Fig. 4.6. For the particular value of F at which the numerical
calculations have been performed, F ≈ 0.0012, the relationship of R versus η appears
to become asymptotic to a power law for larger η. To be more precise the relationship
appears to be
R = 1 + Aη p̃

(4.57)

[where δω is
where the exponent p̃ ≈ 1.26 when F = 0.0012. Moreover, since η = δω

the deviation of the driving frequency from the resonant frequency, see Eq. (2.30)],
we also learn that the relationship R versus δω is asymptotic to a power law at low
δω (i.e. close to resonance) with the exponent −p̃.
Visual inspection of figures such as Fig. 4.3 leads us to notice that an individual
trajectory by itself differs qualitatively very little from a trajectory at F = 0. In other
words, we hypothesize that a cusp-like caustic happens as a result of a nonuniform
perturbation of the zeroth-order trajectories seen in Fig. 4.1, even though the actual
perturbation of each trajectory is small (meaning, for example that δr(t) ¿ r(0) (t)
at a given t). If so, we may be able to capture the phenomenon with a perturbation
theory. Furthermore, it is reasonable to suggest that the perturbation takes place
mostly around r = 1 where trajectories slow down because the symmetry-breaking
component of velocity that appears at nonzero F is most important there,11 as the
nominal velocity is smallest at r = 1. To test this hypothesis we consider a different
system by reversing the sign in the original hamiltonian g [see Eq. (2.40)]: instead

of 41 (Q2 + P 2 − 1) we consider the system with 14 (Q2 + P 2 + 1) . This eliminates
the slowing-down effect. The escape trajectories for this system for η = 0.2 and
F = 0.05 are displayed in Fig. 4.7. We find that caustics can be found in this system
as well. Hence, the slowing-down region around r = 1 is not crucial to the formation
of caustics. In the next subsection we will demonstrate analytically that indeed,
11

In the absence of damping, there is a circle of stagnation at r = 1, while the damping term
introduces a gradient towards the origin.

15
13
11

Radius

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Radius − (Radius at η = 0.03)

−1

Figure 4.6: Radius of the cusp versus η for a fixed F ≈ 0.0012. The asymptotic
exponent p̃ in R = 1 + Aη p̃ is ≈ 1.26.

−2

−4

−6

−8
−8

−6

−4

−2

Figure 4.7: A portion of the caustic in the system based on g(Q, P ) =
(Q2 + P 2 + 1) , indicated by open circles and a sample of 5 trajectories. Note
the absence of the turn-around at r = 1.
the slowing-down region is not responsible, and caustics in fact, arrive from infinity.
It appears that caustics form as a result of the non-uniform perturbation of each
trajectory at large distances.

4.2.4

Analytical Predictions

One way to calculate trajectories is by a straightforward perturbation theory: let
Q(t) = Q(0) (t) + FQ(1) (t), P (t) = P (0) (t) + FP (1) (t), and same for the momenta
pQ and pP , and substitute into the auxiliary EOM, Eqs. (2.51)-(2.54). There will be
four linear equations for first-order corrections Q(1) , etc., but they will have timedependent coefficients. We can partly ameliorate this inconvenience by eliminating
the need for ṗQ and ṗP equations with the help of a separate perturbative analysis of

76
the action S(Q0 , P0 ).12 This will hand us perturbed momenta as a function of Q and
P , which can then be substituted in the remaining ODEs for Q(1) and P (1) .
We would like to have an expression for the action at a point (Q0 , P0 ) as a series
in F, i.e., S(Q0 , P0 ) = S (0) (Q0 , P0 ) + FS (1) (Q0 , P0 ) + ... . Consider a trajectory that
leaves the attractor and arrives to a point (Q0 , P0 ); for simplicity of notation denote
this trajectory by ξ(t). Then the action is given by
Z (Q0 ,P0 )
S(Q0 , P0 ) =

˙ 0 ) dt0 .
L ξ(t0 ), ξ(t

(4.58)

attr

Of course both L and the trajectory itself vary with F, so we should expand both
around F = 0,
"µ
µ (0) ¶
(0)
∂L
∂L
ξ (1) +
L ξ(t), ξ(t)
= L(0) ξ (0) , ξ˙(0) + F
ξ˙(1)
∂ξ (ξ(0) ,ξ̇(0) )
∂ ξ˙ (ξ(0) ,ξ̇(0) )
¡ ¢
+ L(1) ξ (0) , ξ˙(0) + O F 2 .
(4.59)
When integrated, the first term will yield S (0) (Q0 , P0 ) and the term in brackets will
give zero by definition of L(0) : with integration by parts it can be seen as the integral
of the Euler-Lagrange equations based on L(0) along the trajectory which is a solution
to these equations. Hence,
η¡ 2
S(Q0 , P0 ) =
Q0 + P02 + F
= S

(0)

Z t

(Q0 , P0 ) + FS

−∞
(1)

L(1) ξ (0) (t0 ), ξ˙(0) (t0 ) dt0

(Q0 , P0 ).

(4.60)

As long as the zeroth-order trajectories ξ (0) are unique, i.e., they do not reach the
same point at the same time, as is the case in our problem, to calculate S(Q0 , P0 )
we choose that unique zeroth-order trajectory which reaches (Q0 , P0 ), and substitute
for t in the upper limit in Eq. (4.60) the time at which this point is reached. Given
∂S
the scalar field S(Q0 , P0 ) we calculate momenta in the usual way, pQ = ∂Q
, etc.,
12

A better notation is S(Q, P ) emphasizing that S is a scalar field over the (Q, P ) space, but this
may be confused with Q(t) and P (t) of a trajectory, so we will stick to the notation Q0 and P0 , with
the subscript.

77
and we have produced for ourselves momenta to first order in F without the need
to solve any equations! The method can be extended to higher order. We feed the
obtained F-series for momenta into the full equations for Q̇ and Ṗ and proceed with
the conventional perturbation theory, Q(t) = Q(0) (t) + FQ(1) (t), etc., making sure
that all of the terms of a given order are included. At the end, we would obtain
perturbed trajectories (Q(t), P (t)). It may be that nearby trajectories intersect. An
envelope of these intersections is a caustic.
A few words need to be mentioned about the choice of the coordinate system in
which to do the analysis. Polar system seems to be most useful only for obtaining
zeroth-order solutions, because that is when the angular symmetry takes hold. At a
nonzero F, the fixed point moves away from zero and therefore θ(t) is no longer a
small perturbation over its zeroth-order version for large negative times. After some
consideration, it appears that at nonzero F, it is best to perform the analysis in a
Cartesian system in the time domain, i.e. to perform the action method for Q(t) and
P (t) and the associated momenta, rather then r(θ) for example.
For convenience, we restate the zeroth-order solution:
r(0) (t) = r0 eηt ,
θ(0) (t) = θ0 + t −

r02 ¡ 2ηt
e −1 .
2η

Or
Q(0) (t) = eηt [Q0 cos ψ(t) − P0 sin ψ(t)] ,
P (0) (t) = eηt [Q0 sin ψ(t) + P0 cos ψ(t)] ,
Q2 + P02 £ 2ηt
where ψ(t) = t − 0
e −1 .
2η
This solution is designed to arrive to point (r0 , θ0 ) or (Q0 , P0 ) at time t = 0. Remember that the superscript denotes the order of the solution versus time, while the
subscript denotes the point in space; the subscripts were used to avoid confusion with

78
the full solution Q(t) and P (t). Applying the action method we have,
Z 0

Z 0
¡ (0)
(0)
S(r0 , θ0 ) =
r (t), θ (t) dt + F
L(1) r(0) (t), θ(0) (t) dt,(4.61)
−∞
¸ −∞
(ṙ + ηr)2 + rθ̇ − r + r3
(4.62)
L(0) =
1 h³
(1)
= −
rθ̇ − r + r cos θ + (ṙ + ηr) sin θ .
(4.63)
(0)

In view of this, and Eqs. (4.44)-(4.45), see that
S (0) (r0 , θ0 ) =

(1)

η 2
r ,
2 0 Z

(4.64)

r02 ¡ 2ηt
(r0 , θ0 ) = −ηr0
e sin θ0 + t −
e − 1 dt
2η
−∞
µ 2
¶¸¶
Z 1
ξ −1
ln ξ − r0
dξ.
= −r0
sin θ0 +

ηt

(4.65)
(4.66)

Notice that compactly, the dependence of S (1) upon (r0 , θ0 ) appears as parameters
(1)

inside the integral, and pi can be effortlessly obtained by differentiation with respect
to these parameters under the integral sign - we do not have to mess with differentiation of the limits of integrals! We now drop the subscripts such as 0 in r0 or θ0 .
The superscripts indicate the order in the perturbation theory. We have
S = S (0) + FS (1) .

(4.67)

The Jacobian is
µ 2 ¶2
∂ S
∂pQ ∂pP
∂ 2S ∂ 2S
∂pQ ∂pP
J =
∂Q ∂P
∂P ∂Q
∂Q ∂P
∂P ∂Q
· 2 (0) 2 (1)
2 (0) 2 (0)
2 (0) 2 (1)
∂ S ∂ S
∂ S ∂ S
∂ S ∂ S
+F
+ O(β).
∂Q2 ∂P 2
∂Q2 ∂P 2
∂P 2 ∂Q2

(4.68)
(4.69)

The mixed partials of S (0) disappear in view of Eq. (1). Also, in view of this we have
J = η 2 + ηF∇2 S (1) + O(β).

(4.70)

79
The appearance of the Laplacian is fortunate, since it is a geometric object and can
therefore be taken in any coordinates. We will work in polar coordinates, so

∇S

(1)

1 ∂
r ∂r

∂S (1)
∂r

1 ∂ 2 S (1)
r2 ∂θ2

(4.71)

We now compute the derivatives. First note that
r3
∂S (1)
= S (1) +
∂r
So
1 ∂
r ∂r

∂S (1)
∂r

Z 1

(ξ 2 − 1) cos (φ) dξ.

(4.72)

· Z
1 ∂S (1) 1 ∂ r3 1 2
(ξ − 1) cos (φ) dξ .
r ∂r
r ∂r η 0

(4.73)

The first term on the rhs of Eq.(10) can be read off directly from two lines above:
1 ∂S (1)
=−
r ∂r

Z 1

sin (φ) dξ +

Z 1

(ξ 2 − 1) cos (φ) dξ.

(4.74)

(ξ 2 − 1)2 sin (φ) dξ.

(4.75)

The second term is
3r

Z 1

r3
(ξ − 1) cos (φ) dξ + 2

Moreover,
1 ∂ 2 S (1)
r ∂θ

Z 1

Z 1
sin (φ) dξ.

(4.76)

All together,

∇S

(1)

Z 1

r3
(ξ − 1) cos (φ) dξ + 2

Z 1

(ξ 2 − 1)2 sin (φ) dξ.

(4.77)

So the meat of the calculation has to do with understanding the integrals.
Z 1
I1 =

(ξ 2 − 1) cos (φ) dξ,

(4.78)

(ξ 2 − 1)2 sin (φ) dξ.

(4.79)

Z 1
I2 =

80
The analysis of the integrals I1 and I2 , Eqs. (4.78)-(4.79), can be performed with
the use of Γ-functions. Consider

(p)

Z 1

2

ξ −1
ξ exp iθ +
ln ξ − r
dξ
Z 1
iθ iκ
= e e
ξ p ξ α e−iκξ dξ,

(4.80)

, α = ηi and p are either 0, 2 or 4. Next, change variables to x = iκξ 2 .
where κ = 2η

Then

(p)

eiθ eiκ
2(iκ)z

Z iκ

xz−1 e−x dx,

(4.81)

= 2ηi + p+1
. Note that z lies in the right-half z-plane. When r ≫
where z = α+p+1

the upper limit can be replaced by ∞ (notice that z is independent of r, so as r is
increased but η is held fixed, only the upper limit changes, but the structure of the
integrand does not, so its ok to make this replacement). The integral is taken in the
x-plane along the imaginary axis, but the integrand is an analytic function of x, so
we may rotate the contour to go along the real x-axis. Then, the integral coincides
with the definition of the Γ-function in the variable z. Thus,
I (p) =

eiθ eiκ
Γ(z).
2(iκ)z

(4.82)

Notice that in the limit of large η → ∞, z → p+1
, i.e. α → 0. But α = 0 is equivalent
to starting out with the integral that does not contain the ln-term. This justifies
dropping the ln-term in the case of large η.

81
Now,
I (2) I (0)
+ c.c.


i(θ+κ)
2η
2η
 + c.c.,
3/2
1/2
(iκ)
(iκ)
4(iκ) 2η

I1 =

I (4) I (2) I (0)
+ c.c.
2i
i 2i

i(θ+κ)
2η
2η
2η
 + c.c. .
(iκ)5/2
(iκ)3/2
(iκ)1/2
4i(iκ) 2η

(4.83)

I2 =

(4.84)

These approximations match nicely with numerical evaluations of the full I1 and
I2 from Eqs. (4.78)-(4.79). Each term in the expression for I1 is of the form
e 4η ei(θ+κ− 2η − 2 )
4κq

ln κ

qπ

+q ,
2η

(4.85)

. The same is true for each term in the expression for I2 with an extra
where κ = 2η

phase shift by −π/2. Note that both integrals behave like ∼ 1/r for large r. With
the help of Eq. (4.77), the r-dependence of the correction to the momentum (see
Eq. (4.70)) is
F ∇2 S (1) ∼ F r 2 ,

(4.86)

and hence the comparison of this with the zeroth-order part suggests that R ∼ F −1/2 .
The η-dependence is not as simple, but attains simpler forms in the asymptotic
regimes of η ≫ 1 and η ≪ 1. In particular, for large η and r, both I1 and I2
1/2

behave like η r , which implies, with the help of Eqs. (4.77) and (4.70) that
J ∼ η2 + F

r2
η 1/2

(4.87)

82
Equating the two terms in J suggests that
R∼

η 5/4
F 1/2

(4.88)
1/2

For small η and large r, both I1 and I2 behave like η r e 4η Γ

+ 12
2η

Stirling’s asymptotic approxmation of the Γ-function: for large z
Γ(z) ∼
Then

2η

2πz z−1/2 e−z .

. We invoke

(4.89)

2πei(− 2η − 2 + 2η ln 2η ) e− 4η .

(4.90)
1/2

Putting it all together, we see that both I1 and I2 behave like η r ; to our luck,

the exponential factors e− 4η and e 4η cancel each other out and prevent essential
singularity. Thus, we can again conclude that for large r,
R∼

4.3

η 5/4
F 1/2

(4.91)

Physics Without the Slowing-down Region

The only difference is that now the phase has a different sign:
2

ξ −1
ln ξ + r
φ(ξ) = θ +

(4.92)

and therefore

(p)

Z 1

2

ξ −1
ln ξ + r
dξ
ξ exp iθ +
Z 1
iθ iκ
= e e
ξ p ξ αeiκξ dξ,

(4.93)

where κ = 2η
, α = ηi and p are either 0, 2 or 4. Next, change variables to x = −iκξ 2 .

Then

(p)

eiθ e−iκ
2(−iκ)z

Z −iκ

xz−1 e−x dx,

(4.94)

η the upper limit can be

= 2ηi + p+1
. Again, in the limit of r ≫
where z = α+p+1

replaced by −∞ and the contour can be rotated to go to ∞ along the real axis. Thus,
I (p) =

eiθ eiκ
Γ(z),
2(−iκ)z

(4.95)

and the scaling of R with F would remain unchanged.

4.4

Perturbation Theory for r(p)

Our method above was to extract the physics by comparing the first-order correction
to the zeroth order term, as this was a signature when the Jacobian J begins to blow
up. There is no hope of actually observing any singularities of J with the use of the
perturbation theory - for one thing, the first-order correction to the action S (1) is
single-valued, so there is no hope to extract information about the folds directly. On
the other hand, we know from numerical work that for small F the trajectories are
indeed perturbed only a little. The infinity of J is only an artifact of the particular
coordinate system that we use. In reality, the Lagrangian manifold obtains a small
“wrinkle”, which when viewed from just the right perspective has an infinite slope,
but there is nothing special about this perspective. One would suspect then, that
there must exist a perturbation theory which converges and which captures the real
perturbed trajectories.
One such possibility is to express the manifold as r(p) rather then p(r), which
means we have to express the action as a function of momenta, call it S̃(p). Since
p = ∇S(r), S̃(p) is just the Legendre transform of S(r). In other words
S̃(pP , pQ ) = pP P (pQ , pP ) + pQ Q(pQ , pP ) − S(Q(pP , pQ ), P (pQ, pP )).

(4.96)

84
This is true, but we would like an expression for S̃(pP , pQ ) of the form of Eq. (4.65).
We follow a procedure similar to the procedure employed in Landau and Lifshitz to
perform the Legendre transformation on the Lagrangian and to derive the Hamiltonian. Illustrated on a 1-variable action we have:
∂S
∂S
dq +
dt
∂q
∂t
= pdq − Hdt.

dS =

(4.97)

Add qdp to both sides of the equation. Then
dS + qdp = d(qp) − Hdt.

(4.98)

d(S − qp) = −(qdp + H)dt,

(4.99)

Rearranging terms we have

and we recognize the lhs as the differential of S̃. Hence
S̃(p) = −

q dp −

H(q(t), p(t)) dt

(4.100)

= −

q(t)ṗ(r) + H(q(t), p(t)) dt.

(4.101)

Now, let the perturbed trajectory be given by r(t) = r(0) (t) + Fr(1) (t) and similar for
momenta p(t) = p(0) (t) + Fp(1) (t). They have the same initial and final points in the
phase space as the unperturbed trajectories. Then
S̃ = −

Z nh
ih
q (0) + Fq (1) ṗ(0) + F ṗ(1) dt

´
(0) ¯
(0) ¯
∂H
∂H
dt.
q (1) (t) + F
p(1) (t) + FH(1) q (0) , p(0)
+ H(0) (q (0) , p(0) ) + F
∂q ¯ (0) (0)
∂p ¯ (0) (0)

(4.102)

85
With the help of the canonical equations we have
S̃ = S

(0)

−F

£ (1) (0)
¢¤
q ṗ + q (0) ṗ(1) − q (1) ṗ(0) + q̇ (0) p(1) + H(1) q (0) , p(0) dt. (4.103)

The first and third terms obviously cancel out and we can perform integration by
parts on the fourth term. The boundary terms are zero because p(1) is zero at the
initial and final times and the other term cancels out the second term in the integrand
above. We are this left with
S̃ = S

(0)

−F

H(1) q (0) , p(0) dt,

(4.104)

which now needs to be reexpressed in terms of momenta. Recall that at zeroth-order
in F we have pQ = ηQ and pP = ηP . So we can just use Eq. (4.64) to get
S̃ (0) (pQ , pP ) =

1 ¡ 2
pQ + p2P .
2η

(4.105)

Now,
H(1) = pP .
Recall also that

(0)

r2 ¡ 2ηt
= re sin θ + −
e −1 ,
2η

(4.106)

ηt

(4.107)

(0)

and as was just mentioned, pP = ηP (0) . Therefore, the first-order term in Eq. (4.104)
is exactly the same as that given by Eq. (4.65), except that the result must be expressed in terms of pP and pQ instead of P and Q, but that is trivial with the relations
just given. Eq. (4.99) tells us that q = ∇p S̃ - this is the alternative expression for

86
the manifold q(p). The new Jacobian is
2 (0) 2 (1)
2 (1) 2 (0)
2 (0) 2 (0)
+F
J˜ =
∂p2Q ∂p2P
∂p2Q ∂p2P
∂p2Q ∂p2P
F ∂S (1) ∂S (2)
= 2+ 3
∂Q2
∂P 2
= 4.
Thus, the new method does not predict anything new!

Part II
Renormalization Group Method
for Predicting Frequency Clusters
in a Chain of Nearest-Neighbor
Phase Oscillators

Chapter 5
Synchronization: Background
More then three centuries have passed since Huygens’s observation of entrainment of
clocks [90]. In modern days, the study of synchronization in nature has been a subject
of major revival, yet some of the most challenging questions remain unanswered.
A variety of models have been in use for these systems. A “toy model” for exploration of various issues in this field is the phase model, which represents oscillators
that interact with each other through an odd function of their phase differences,
usually taken to be sine:
θ̇i = ωi +

Ki,j sin (θj − θi ).

(5.1)

i6=j

Not unlike the Ising model of magnetism, this phase model of oscillators is sufficiently
simple to allow for the development of intuition, yet it contains enough physics to
capture some of the phenomenology exhibited by more realistic models. An analysis
of a more general model that allows for a change of an amplitude of each oscillator
reveals that under certain conditions the dynamics develops an attracting limit cycle
with a fixed (or approximately fixed) amplitude, and can thus be reduced to the
aforementioned phase model [47, 91]. An example of a procedure like this has been
given recently in the theoretical study of nanomechanical oscillators [92]. In the phase
mode the set of frequencies {ωi } is random and is sampled from a distribution g(ω).
Analogously to systems in equilibrium Statistical Mechanics, questions about the
transition to global synchronization arises in the context of synchronization as well; in

89
the limit N → ∞ one may expect abrupt transitions to globally synchronized states.
As in the equilibrium realm, such a transition is reflected in the abrupt change in an
order parameter. For oscillators, one choice of an order parameter is

1 X iθj
e ,
N j=1

(5.2)

which may attain values À 1/ N when some macroscopic fraction of oscillators
shares a common period. The question of the existence of such a transition to macroscopic synchronization is subject to dimensionality, properties of coupling Ki,j and
properties of g(ω). For example, the well-known all-to-all coupling Ki,j = K/N
treated by Kuramoto [46],
θ̇i = ωi +

KX
sin (θj − θi ),
N i6=j

(5.3)

allows for a threshold value of K = Kc beyond which an order parameter attains
values À 1/ N . Due to the mean-field nature of the all-to-all interaction, the value
of Kc can be found analytically - it is 2/πg(0) [45, 47, 93]. Other features of this
transition to the synchronized state, such as the growth of the fraction of locked
oscillators as the coupling increases above the critical value, can also be calculated.
Many systems where synchronization might occur involve short-range coupling
[48, 49], such as oscillators on a lattice with nearest-neighbor interactions. Although
nanomechanical examples have not been constructed yet, it is plausible that interactions will be short range if carried by the support or a substrate between oscillating
nanobeams. Systems with short range coupling are much harder to analyze, and much
of our understanding of these systems rests on numerical simulations of the equations
of motion as there are few theoretical tools that yield a quantitative description of
the collective motion.
The focus here will be on a one-dimensional chain of phase oscillators interacting

90
with their nearest neighbors.
θ̇i = ωi + K sin (θi−1 − θi ) + K sin (θi+1 − θi ).

(5.4)

In spite of its simplicity, this model presents a challenge as we will see below. It
has been shown by Strogatz and Mirollo that a nearest-neighbor chain cannot exhibit extensive synchronization for any finite vale of K [51]. Strogatz and Mirollo
also described the probability for a global synchronization at finite N and how this
probability tends to zero as N → ∞ [52]. Although global synchronization is not
possible, interesting collective structures do occur. In particular, oscillators tend to
organize themselves into clusters of common frequency, ω̄, defined as
ω̄i ≡

θi (t) − θi (t0 )
(t−t0 )→∞
t − t0
lim

(5.5)

such that each oscillator in a given cluster shares its ω̄ with all the other oscillators in
the same cluster [51], but the dynamics of different oscillators in a given cluster may
be different. Hence, the structure of the unsynchronized phase in the nearest-neighbor
chain is quite complex. The relevant questions are regarding the distribution of cluster
lengths and distribution of cluster frequencies as a function of the distribution g(ω)
of frequencies.1 {ωi } As the coupling strength K is increased, the average size of such
clusters grows, Fig. 5.1, but it does not become infinite at a finite K [51].
In a sense, understanding the properties of such clusters is a more complicated
problem then predicting the feasibility of global synchronization, because the latter
question is equivalent to predicting the existence of an attracting fixed point of the
dynamics (this by itself is far from trivial in the N → ∞ limit), while there is no such
clear strategy in the case of frequency clusters - they are a dynamic phenomenon.
Twenty years ago, Strogatz and Mirollo wrote [51]:
“The dynamical behavior of θ̇i = ωi +K

hi,ji sin (θj − θi ), j ∈ {1, ..., N }

is not well understood in the regime before phase-locking occurs. In par1

well

also as a function of the distribution G(K) of couplings {Ki } if we allow disorder in these as

Figure 5.1: Development of frequency clusters as the value of coupling K is increased,
shown on a small section of a 1000 oscillator chain. Solid diamonds: time-averaged
running frequencies of oscillators, defined similarly to Eq. (6.40) with t − t0 a large,
but finite time interval; lines: intrinsic frequencies that appear in Eqn. (5.4).
ticular, it is not known if or how the distribution of number and size of
synchronized clusters scale with K, N , and d”
where d is the dimensionality of the system.
The understanding of these structures from the point of view of dynamical systems
has been deepended by a series of papers by by Ermentrout and Kopell in early 1980s.
For example, in [94] they considered a chain of oscillators with frequencies laid out
linearly (not randomly) along the chain, proved the existence of limit cycles in the
θ-space
and related the first cluster break in the chain as the coupling strength is
lowered to the appearance of a limit. Using this they were able to predict sizes of
frequency differences between neighboring clusters. Around the same time, in a rather
unrelated field, Dasgupta and Ma [95, 96], developed a renormalization group (RG)
analysis of the ground state of 1-dimensional random quantum spin chains (see also
more recent work by Daniel Fisher and co-workers [97, 98]). Their scheme was based
on solving for the energy eigenvalues of the small parts of the chain (two spins), while
treating the effect of the neighbors via a secon-order perturbation theory. Using this,

92
they were able to progressively thin out the spins interacting via largest couplings
and in this way to set up an RG flow.
The main idea behind the strong disorder approach applied to quantum-mechanical
problem is best demonstrated in the case of the spin-1/2 Heisenberg model:
Ĥ =

Ji Ŝi · Ŝi+1 .

(5.6)

The ith term in the sum considered on its own, splits the four states of the spins i
and i + 1, into a singlet state, and three triplet states excited by the energy Ji . Thus
a real-space decimation step can be proposed: freeze the spin pair with the strongest
interaction, Jn = max{Ji }, into a singlet state. Perturbative corrections to this state
introduce an effective coupling between the neighbors of the spins n and n + 1 which
is:
HRG =

Jn−1 Jn+1
Ŝn · Ŝn+1 ,
2Jn

(5.7)

an interaction identical in form to the operators in the bare Hamiltonian, but with
Jn+1
< Jn−1 , Jn , Jn+1 . By repeating this
a much suppressed strength: Jef f = Jn−1
2Jn

decimation step the number of free spins is gradually decreased; at the stage when all
spins are bound into singlets, we obtain the ground state of the Hamiltonian. Note
that this procedure is possible if the coupling to the rest of the chain is much weaker
then between the two spins in question. Thus, their method was best applicable to
chains that contained a large fraction of couplings which are significantly different
from adjacent couplings, in other words, this technique applied to chains with broad
coupling distributions.
More recently, the technique of strong randomness RG has been extended to explain the superfluid-insulator transition in random bosonic chains [99, 100] (a quantum reactive analog of the nearest-neighbor chain, where the θ̈ replaces θ̇, and θ is a
quantum operator). This latter system consists of superconducting grains of random
sizes coupled to each other via Josephson couplings which are also random. Those
superconducting grains with exceptionally large charging energies tend to introduce a

93
break in the chain, but perturbative calculation demonstrates that the effect of such
a grain is to establish a weak effective coupling between its adjacent grains. Those
couplings which are exceptionally large tend to lock together the phases of the order
parameter of the superconducting grains which they couple, thus allowing these two
grains to be treated as one grain with a new effective charging energy. An RG based
on these decimation steps was built and its analysis led to the mapping out of a
phase diagram for parameters under which the system behaves as a superconductor
and those under which it behaves as a Mott insulator. This served as our motivation
for an analogous task with the chain of oscillators.
The rest of the work serves a dual purpose - to develop an RG method as a quick
and efficient tool for the analysis of the nearest-neighbor model and then use this tool
to study the physics of frequency clusters exhibited by this model. The content of
this part of the thesis follows the recent paper by Kogan, Rogers, Refael, and Cross
[101], although the presentation is significantly expanded. The presentation in [101]
is less intensive on the algebra, and therefore forms a good supplement to the more
detailed derivations here, especially in Sections 6.2 and 6.3.

Chapter 6
Oscillator Chain RG
6.1

Idea

We will develop the strong-disorder real space RG method to analyze the random oscillator chain. The nearest-neighbor chain presents a classical problem with dissipative
equations of motion, rather than a quantum problem with a conserved Hamiltonian.
Nevertheless, we can apply the strong disorder RG method to a classical problem as
well:
θ̇i = ωi + Ki−1 sin (θi−1 − θi ) + Ki sin (θi+1 − θi ),

(6.1)

where we generalized to introduce the randomness into the couplings. The distributions of ωi is given by g(ω) and couplings Ki come from the distribution G(K).
These random frequencies and couplings are treated roughly speaking, as energies in
a Hamiltonian. Indeed, a strong coupling between two oscillators will tend to force
them to rotate at the same frequency, which allows for a renormalized description of
two oscillators by one effective oscillator. Concomitantly, oscillators with frequencies
which are anomalously different from each of their neighbors’ frequencies will rotate
essentially independently, or as free oscillators, with only small perturbative corrections from their neighbors; these neighbors are in turn perturbed only slightly by the
anomalous oscillator. This allows us to get rid of such anomalous oscillators from
the description of the chain and replace it with the description in which each of the
neighbors’ parameters are slightly changed.

95
As the RG proceeds the parameters Ki and ωi are modified and the strength
of their disorder changes.1 As equations of motion are renormalized, they generally
acquire new terms. We list these terms, but ignore them in this work, except for the
appearance of the parameter mi :
mi θ̇i = mi ωi + Ki−1 sin (θi−1 − θi ) + Ki sin (θi+1 − θi ).

(6.2)

It will follow from the strong-coupling decimation step that when two oscillators
with mi and mi+1 are described by one effective oscillator, its m is given simply by
mi +mi+1 . Hence, the value of m equals the number of original oscillators represented
by this effective oscillator. The original chain, which is the starting point for RG will
have all mi = 1. It will become apparent that oscillators with larger m are less
affected by perturbations from their neighbors, hence m is a measure of sensitivity
of an oscillator to outside perturbations (similarly to inertia in Mechanics) and will
be referred to as its “mass”, although it is not a true mass in the sense of Newton’s
second law. We now proceed to develop each decimation step in detail.

6.2

Strong-Coupling Step

Figure 6.1: A small piece of the chain around the strongest coupling, labelled here as
Kn . All the neighboring Ki and all the local ωi are assumed to have much smaller
magnitude then Kn .
Consider the largest coupling in the chain: suppose this coupling is Kn . Because of
strong disorder, the frequencies ωn and ωn+1 of the pair of oscillators that Kn couples
as well as the couplings Kn−1 and Kn+1 to the neighbors of the pair will almost

Even if the bare chain contains all identical Ki , the RG will create also Ki = 0. However, we
do not expect the RG to work well on chains with all identical Ki .

96
always be small compared to Kn . Associated with this, there are small parameters
ωn /Kn , Kn−1 /Kn , etc. We will occasionally use a symbol ² to denote quantities of
this order. Intuitively we know that oscillators n and n + 1 rotate together as one
oscillator, with a small phase difference between them that remains bounded. Notice
that it is not the large magnitude of Kn by itself which causes this; the phase difference
between θn and θn+1 can not be assumed to be bounded if the coupling Kn between
them is very large but the difference between their intrinsic frequencies or another
neighboring Kn are comparably large. Therefore, it is the small value of ratios such
as ωn /Kn and Kn−1 /Kn that determine whether an oscillator is subjected to the
strong-coupling step.
Because we expect the phase difference to be small, it will be convenient to switch
from variables θn and θn+1 to the mass-weighted average phase Θ and the phase
difference variable δ defined as
Θ ≡ (mn θn + mn+1 θn+1 )/M,

(6.3)

δ ≡ θn+1 − θn ,

(6.4)

and

respectively, where M = mn + mn+1 . For reference, the equations of motion for the
strongly coupled oscillators and their neighbors read:
mn−1 θ̇n−1 = mn−1 ωn−1 + Kn−2 sin (θn−2 − θn−1 ) + Kn−1 sin (θn − θn−1 ), (6.5)
mn θ̇n = mn ωn + Kn−1 sin (θn−1 − θn ) + Kn sin (θn+1 − θn ),

(6.6)

mn+1 θ̇n+1 = mn+1 ωn+1 + Kn sin (θn − θn+1 ) + Kn+1 sin (θn+2 − θn+1 ),

(6.7)

mn+2 θ̇n+2 = mn+2 ωn+2 + Kn+1 sin (θn+1 − θn+2 ) + Kn+2 sin (θn+3 − θn+2 ).(6.8)
Rewriting equations for θn and θn+1 in terms of Θ and δ we get
mn+1 ´
mn ´
M Θ̇ = M ω 0 + Kn−1 sin θn−1 − Θ +
δ + Kn+1 sin θn+2 − Θ −
δ . (6.9)

97
Similarly, oscillator θn+2 is coupled to Θ and δ

mn ´
mn+2 θ̇n+2 = mn+2 ωn+2 + Kn+1 sin Θ − θn+2 +
δ + Kn+2 sin (θn+3 − θn+2 ),
(6.10)
and analogously for oscillator θn−1 . The phase difference δ satisfies
δ̇ = (ωn+1 −ωn )−

mn ´ Kn−1
mn+1 ´
Kn
Kn+1
sin θn+2 − Θ −
sin θn−1 − Θ +
sin (δ)+
δ −
δ .
mn+1
mn
(6.11)

By the strong-disorder assumption, the term sin δ has a much larger prefactor Kµn then
the other terms, so this term is responsible for quick relaxation of δ to the quasi-static
equilibrium value given by
¶2
Kn−1
µ Kn+1
sin (θn+2 − Θ) −
sin (θn−1 − Θ) + O
, (6.12)
δ(t) = δ0 +
Kn mn+1
mn
Kn
−ωn
where δ0 = µ ωn+1
. So δ consists of small oscillations of order ² centered around
Kn

the time-independent mean δ0 which is also of order ² and is proportional to ωn+1 −
ωn . The phase δ will mediate an interaction between oscillators n − 1 and n + 2,
since each of these oscillators is coupled to δ. When this solution is substituted into
n+1
δ we get
Kn−1 sin θn−1 − Θ + mM
mn+1 ´
Kn−1 sin θn−1 − Θ +
δ0
µ Kn−1 Kn+1
mn+1 Kn−1
(sin (θn+2 + θn−1 − 2Θ) + sin (θn+2 − θn−1 )) −
sin (2θn−1 − 2Θ)
2Kn
mn 2Kn
(6.13)

and when it is substituted into Kn+1 sin θn+2 − Θ − mMn δ we get
mn ´
Kn+1 sin θn+2 − Θ −
δ0
µ Kn−1 Kn+1
mn Kn+1
(sin (θn+2 + θn−1 − 2Θ) + sin (θn−1 − θn+2 )) −
sin (2θn+2 − 2Θ) .
2Kn
mn+1 2Kn
(6.14)

98
Therefore, Eq. (6.9) transforms into
mn+1 ´
mn ´
δ0 + Kn+1 sin θn+2 − Θ −
δ0
M Θ̇ = M ω 0 + Kn−1 sin θn−1 − Θ +
µ Kn−1 Kn+1
mn+1 Kn−1
sin (θn+2 − θn−1 − 2Θ) −
sin (2θn−1 − 2Θ)
Kn
mn 2Kn
mn Kn+1
sin (2θn+2 − 2Θ) , (6.15)
mn+1 2Kn
where Ω is the mass-weighted frequency,
ω0 =

ωn mn + ωn+1 mn+1

(6.16)

M = m1 + m2 ,

(6.17)

and

while Eq. (6.10) transforms into
mn ´
mn+2 θ̇n+2 = mn+2 ωn+2 + Kn+1 sin Θ − θn+2 +
δ0 + Kn+2 sin (θn+3 − θn+2 )
mn Kn+1
µ Kn−1 Kn+1
(sin (2Θ − θn+2 − θn−1 ) − sin (θn−1 − θn+2 )) −
sin (2Θ − 2θn+2 ) .
2Kn
mn+1 2Kn
(6.18)

An analogous expression exists for θn−1 . To leading order it turns out that the oscillating corrections to δ give small extra coupling terms of forms not included in the
equation of motion Eq. (6.2): a second-nearest-neighbor interaction between oscillators, a three body interaction involving the new effective oscillator and the neighboring pair, and a second harmonic correction to the form of the pairwise interaction.
As is conventional in real space renormalization group approaches, we will neglect
these more complicated interaction terms. Then, aside from the factor δ0 resulting
equations take the original form. Fortunately, in one dimension the phase δ0 can be
unwound by redefining the phases θn−1 and all the phases prior to it as
θj0 = θj +

mn+1
δ0 , j ≤ n − 1,

(6.19)

99
as well as the phase θn+2 and all the phases following it as
θj0 = θj −

mn
δ0 , j ≥ n + 2.

(6.20)

Finally, we arrive at
mn−1 θ̇n−1 = mn−1 ωn−1 + Kn−2 sin (θn−2 − θn−1 ) + Kn−1 sin (Θ − θn−1 ),(6.21)
M Θ̇n = M ω 0 + Kn−1 sin (θn−1 − Θ) + Kn sin (θn+2 − Θ),

(6.22)

mn+2 θ̇n+2 = mn+2 ωn+2 + Kn+1 sin (Θ − θn+2 ) + Kn+2 sin (θn+3 − θn+2 ).(6.23)
In summary, the effect of the strong-coupling decimation step is to replace the pair
of oscillators with phases θn and θn+1 by a single oscillator with phase Θ, frequency
ω 0 and mass M respectively given by Eqs. (6.3), (6.16), and (6.17).

6.3

Crazy Oscillator Step

Figure 6.2: A small piece of the chain around the oscillator with the largest frequency,
labelled here as Ω. All the neighboring Ki and all the neighboring ωi are assumed to
have much smaller magnitude then ωn .
Consider the oscillator with the largest frequency in the entire chain: suppose this

100
is oscillator n. We use Ω for its frequency ωn to emphasize this is a large quantity.
Because of strong disorder the coupling strengths Kn−1 and Kn to the neighboring
oscillators, and the frequencies ωn±1 of these oscillators will almost always be small
compared to Ω. Associated with this, there are small parameters Kn−1 /Ω, ωn−1 /Ω,
etc. In the present section we will use the symbol ² to denote quantities of this
order. To zeroth order in ² the nth oscillator runs freely at the frequency Ω. We call
such oscillators “crazy oscillators”, Fig. 6.2. However, as in the case of the strong
coupling, it is not the absolute magnitude of Ω that determines whether an oscillator
is subjected to this step, but the O(²) ratios such as Kn−1 /Ω and ωn−1 /Ω: the phase
of an oscillator with a large Ω, but strongly coupled to a slow neighbor or positioned
next to another fast neighbor can not be assumed to advance freely (see also Section
6.3.1 for elaboration).
Deviations from the free-running solution can be calculated perturbatively in ².
The most important effect we might look for is an induced effective coupling via
oscillator n between the neighbors n − 1 and n + 1 (similar to what we have seen in
the strong-coupling step), since this might lead to their synchronization. However,
at least to order ²2 as we describe below, we find that there is no induced interaction
of the type that tends to lock the neighboring oscillators. 2 Thus the main effect of
decimating the largest frequency oscillator in the chain is to eliminate this oscillator
from further consideration (it forms a cluster in the final tally at frequency close to Ω
and with a size given by the mass parameter m) and to cut the chain into two pieces.
This is a key feature of the strongly random chain that limits the size of synchronized
clusters. At order ²2 we do find a renormalization of the frequencies of both the
eliminated oscillator and the two neighbors remaining in the chain. We now analyze
the elimination of the crazy oscillator (oscillator n in Fig. 6.2) perturbatively. The

This should be compared with the strongly disordered spin chain, where an effective interaction
between the neighbors is given by eliminating the high energy spins.

101
EOM for oscillator n and the neighboring oscillators read:
dθn−1
Kn−1
Kn−2
sin (θn − θn−1 ) +
sin (θn−2 − θn−1 )
= ωn−1 +
dt
mn−1
mn−1
dθn
Kn−1
Kn
= Ω+
sin (θn−1 − θn ) +
sin (θn+1 − θn )
dt
mn
mn
Kn+2
dθn+1
Kn+1
sin (θn − θn+1 ) +
sin (θn+2 − θn+1 )
= ωn+1 +
dt
mn+1
mn+1

(6.24)

(6.25)

(6.26)

As Ω/Kn and Ω/Kn−1 become very large, the dynamic of a crazy oscillator is given
approximately by θn = Ωt. The corrections are caused by the interaction with the
neighbors. There will be corrections to the average frequency at which θn grows,
but there will also be new periodic terms that arise due to the periodic slow down
and speeding up as the phase θn advances in the “washboard potential” of the slow
neighbors. There will be other effects, for example, the neighbors pick up a fast
component to their dynamics imparted by the oscillator n and this change will be
reflected back in the dynamics of the nth oscillator itself. Regardless of the effect, we
must make the ansatz that allows for both the change in the average growth rate of
θn and the appearance of small new (periodic) terms. Thus, we let
= Ω 1 + + 2 + ... t + + 2 + ...
Ω Ω
Ω Ω
= Ω0 t + + 2 + ... .
Ω Ω

θn

(6.27)

This is the perturbation ansatz for θn . Occasionally we will write Ω0 to save space
rather then the series for which it stands. As for the dynamics of each of the neighbors,
θn±1 , let l(t) be the solution to n − 1st oscillator when it is completely disconnected

102
from the crazy oscillator, i.e., this l(t) satisfies:
dl
Kn−2
sin (θn−2 − l)
= ωn−1 +
dt
mn−1
dθn−2
Kn−3
Kn−2
= ωn−2 +
sin (θn−3 − θn−2 ) +
sin ((l − θn−2 )
dt
mn−2
mn−2
...

(6.28)
(6.29)

(an analogue of l(t) for n + 1st oscillator will be called r(t)). Then, at finite Kn−1 /Ω,
the dynamic of oscillator n − 1 is given by
Kn−1
θn−1 = l(t) −
cos (Ωt − l(t)) + O
mn−1 Ω

Ω2

(6.30)

as can be verified by direct differentiation and substitution of the crazy oscillator
ansatz, Eq. (6.27) and the definition of l(t), Eq. (6.30). This gives the fast correction
imparted onto the slow motion of a nearest-neighbor of a crazy oscillator, and we will
calculate the effect of this fast component of the motion of the neighbor back on CO.
We will now take the self-consistent, iterative approach to extract the first few
terms in the perturbation series for θn . We have
´ ψ̇
= Ω 1 + + ... + + ...
´ ´ K
´ ´
Kn−1
= Ω+
sin l(t) − Ω 1 + + ... t +
sin r(t) − Ω 1 + + ... t .
mn
mn
(6.31)

θ̇n

We made an assumption that ψ/Ω, etc. terms are much smaller then the terms we
keep in the arguments of sines (we can not assume Ωα + ... t to be small because
t can be large). We will calculate ψ and verify that the “smallness” assumption
was correct, hence making the method self-consistent. Then integrating the equation
while treating l(t) and r(t) as constants we have

µ ¶
Kn−1
Kn
+ 2 + ... t + =
cos l(t) − Ω t +
cos r(t) − Ω t + O
. (6.32)
Ω Ω
mn Ω
mn Ω
Ω2

Comparing terms we see that α must be set to zero, since it was defined to be a

103
constant, not a function of time. Then we have found the first order corrections ψ:
¶ ¶
¶ ¶
Kn−1
Kn
cos l(t) − Ω 1 + 2 + ... t +
cos r(t) − Ω 1 + 2 + ... t .
ψ=
mn
mn
(6.33)
We now repeat the procedure:
ψ̇
φ̇
Ω 1 + 2 + ... + + 2 ... =
Ω Ω
Kn−1
Kn−1
Kn
Kn−1
sin l(t) −
cos (l(t) − Ω0 t) − Ω0 t −
cos (l(t) − Ω0 t) −
cos (r(t) − Ω0 t)
Ω+
mn
mn−1 Ω
mn Ω
mn Ω
Kn
Kn
n−1
sin r(t) −
cos (r(t) − Ω0 t) − Ω0 t −
cos (l(t) − Ω0 t) −
cos (r(t) − Ω0 t)
mn Ω
mn+1 Ω
mn Ω
mn Ω
+ ... .

(6.34)

Again, in the spirit of the self-consistent approach, we neglect φ-terms from the
arguments of the sines. We now expand the sines in 1/Ω and cancel out the zerothorder terms on the right-hand side with the ψ̇ term on the left-hand side, leaving
φ̇
+ ... + 2 + ...
Kn−1
Kn−1 Kn
=−
(1 + cos (2l(t) − 2Ω0 t)) −
cos (l(t) − Ω0 t) cos (r(t) − Ω0 t)
2mn µn,n−1 Ω
m2n Ω
Kn Kn−1
Kn2
(1 + cos (2r(t) − 2Ω0 t)) −
cos (r(t) − Ω0 t) cos (l(t) − Ω0 t) + ...
2mn µn,n+1 Ω
m2n Ω
Kn−1
Kn2
(1 + cos (2l(t) − 2Ω0 t)) −
(1 + cos (2r(t) − 2Ω0 t))
2mn µn,n−1 Ω
2mn µn,n+1 Ω
Kn Kn−1
[cos (l(t) + r(t) − 2Ω0 t) + cos (r(t) − l(t))] + ...,
(6.35)
m2n Ω

=−

−1
−1
with µ−1
i,j = mi + mj . From this we learn two things. First, we read off β:

β=−

Kn−1
Kn2
2mn µn,n−1 2mn µn,n+1

(6.36)

104
Second, we read off an equation for φ:
Kn−1
φ̇
Kn2
cos (2l(t) − 2Ω0 t)) −
cos (2r(t) − 2Ω0 t)
= −
2mn µn,n−1
2mn µn,n+1
Kn Kn−1
Kn Kn−1
cos (l(t) + r(t) − 2Ω0 t) −
cos (r(t) − l(t))
(6.37)
mn
m2n

We will not proceed with integrating this; what we learned at this stage is that an
equation for the crazy oscillator can be rewritten without explicitly containing the
variable θn , at the expense of introducing extra terms:

Kn−1
Kn2
Kn Kn−1
θ̇n = Ω 1 −
cos (r(t) − l(t))
2mn µn,n−1 Ω2 2mn µn,n+1 Ω2
mn2 Ω2
¢ Kn
Kn−1
sin l(t) − Ω0 t +
sin r(t) − Ω0 )t
mn
mn
Kn−1
Kn2
cos 2l(t) − 2Ω0 t ) −
cos 2r(t) − 2Ω0 t
2mn µn,n−1 Ω
2mn µn,n+1 Ω
Kn Kn−1
cos (l(t) + r(t) − 2Ω0 t) + ...,
m2n Ω

where

Ω = Ω 1 + 2 + ... ,

(6.38)

(6.39)

and β is given by Eq. (6.36).
We can define a long-time average rate at which θn advances as
ω̄n ≡ hθ̇n it =

Z T
θ̇n dt =
t0

θi (t0 + T ) − θi (t0 )

(6.40)

This is very similar to the definition in Eq. (6.40), with one difference that T is some
characteristic time of the slow dynamics (for example of r(t) and l(t)). Over this time
period, the parts that depend explicitly on Ω0 t are much suppressed and we have
Kn−1
Kn2
Kn Kn−1
cos (θn+1 − θn−1 ) ,
ω̄n ≈ Ω 1 −
2mn µn,n−1 Ω2 2mn µn,n+1 Ω2
m2n Ω2

(6.41)

where we have also substituted l(t) → θn−1 and r(t) → θn+1 without introducing
errors at the order considered. This is our final result for the crazy oscillator: we

105
started with the unperturbed motion given by θ̇n = Ω and we calculated the longtime average after taking the effect of the neighbors into account to second order in
1/Ω, resulting in Eq. (6.41). Stepping back we see that there are two O(²2 ) effects
that slightly perturb the average rate of the uniform advancement of the phase θn .
The first is the fact that crazy oscillator rotates in the tilted washboard potential
created by its neighbors. Since typically ωn−1 , ωn+1 ¿ Ω for strong disorder, this
potential can be treated as effectively static. The second effect is a small component
of fast dynamics added to the motion of the neighboring oscillators, which then acts
back on the crazy oscillator. Both of these effects change the average frequency of the
crazy oscillator by a shift of order K 2 /Ω. The term cos (r(t) − l(t)) comes only from
the first effect - the positions of each of the neighbors affects washboard potential,
and hence, the slowing down and speeding up of the crazy oscillator.3
Corresponding corrections happen to the motion of the neighboring oscillators
as well. Indeed, a similar procedure can be performed for either of CO’s neighbors
and we can obtain a renormalized equation for both θn±1 which will not be explicitly
coupled to θn . Thus, we will remove θn from the description of the chain! For example,
when we substitute perturbation anzatz for θn , Eq. (6.27) and for θn−1 , Eq. (6.30)
(and analogously θn+1 ), then Eq. (6.24) for the dynamics of θn−1 becomes, after all
the algebra
Kn−1
Kn−1 Kn
cos (l(t) − r(t))
2mn µn,n−1 Ω 2mn mn−1 Ω
Kn−1
Kn−1
sin (Ω t − l(t)) −
cos (2l(t) − 2Ω0 t)
mn−1
2mn−1 µn,n−1 Ω
Kn Kn−1
Kn−2
cos (l(t) + r(t) − 2Ω0 t) +
sin (θn−2 − θn−1 ) + ... .
2mn−1 mn Ω
mn−1
(6.42)

θ̇n−1 = ωn−1 +

As in the case of the crazy oscillator n, we drop the terms containing Ω0 t because the
effect of these terms is to introduce small oscillations to the dynamics θn−1 on top of

The somewhat different derivation which is based directly on calculating these two aforementioned effects, rather then a more formal manipulation of the equations of motion is presented in
[101].

106
whatever dynamics this oscillator has without these terms (of course oscillator θn−1
is coupled to other oscillators, for example to θn−2 , so this nominal dynamics may be
complicated). We then again make the substitution l(t) → θn−1 and r(t) → θn+1 to
end up with
θ̇n−1 = ωn−1 +

Kn−1
Kn−1 Kn
Kn−2
cos (θn−1 − θn+1 )+
sin (θn−2 − θn−1 )+... .
2mn µn,n−1 Ω 2mn mn−1 Ω
mn−1

(6.43)

Finally, we argue that the cosine term that couples oscillators n − 1 and n + 1 can
be dropped. It turns out that the corresponding term in the equation for θ̇n+1 has
same sign! In other words, this new type of interaction does not have action-reaction
symmetry. This sign peculiarity arises because the decimated oscillator has a particular direction of rotation. Due of the lack of action-reaction symmetry, the cosine
term does not tend to pull neighboring oscillators together, and so does not directly
influence the clustering. We therefore omit this term in the renormalization procedure. We will also ignore this term in the motion of the crazy oscillator Eq. (6.41)
since it averages to zero on the long timescale of the period of oscillators n ± 1 [longer
timescale then the characteristic time T over which the average was performed to
arrive at Eq. (6.41)].
So far, the crazy oscillator decimation step can be summarized as follows: (i) remove the crazy oscillator from the chain and set the interaction between its neighbors
to zero; (ii) adjust the frequencies of each of the neighbors by
δωn−1 =

Kn−1
Kn2
, δωn+1 =
2mn−1 µn,n−1 ωn
2mn+1 µn,n+1 ωn

(6.44)

with µ the reduced mass of the pair; (iii) adjust the frequency of the crazy oscillator

107
itself by the corresponding correction from each of the neighbors4
Kn−1
Kn2
δ ω̄n = −
2mµn,n−1 Ω 2mn µn,n+1 Ω

(6.45)

Two observations are in order. First, notice that the constant frequency renormalizations of the crazy oscillator is reciprocal to the sum of its neighbors. This is because
at the pairwise level interactions between the oscillators conserve the mass-weighted
frequency average: by adding all equations of motion we can see that N
i=1 mi δωi = 0,
which is just a consequence of the fact that interactions are odd. Moreover, at the
level of approximation that we consider, δωi = 0 for all but the crazy oscillator n and
its nearest neighbors, n − 1 and n + 2. Hence mn−1 δωn−1 + mn δωn + mn+1 δωn+1 = 0.
Second, the fact that frequency correction terms (the first two terms in Eq. (6.41)
are additive suggests that they can be obtained from the analysis of two 2-oscillator
systems, each of which can be solved exactly. The phase difference φ = θn − θn−1 of
two coupled oscillators satisfies
φ̇ = (ωn − ωn−1 ) −

Kn−1
sin φ.
µn,n−1

(6.46)

n−1
When ωn − ωn−1 ≤ µKn,n−1
the phase difference φ is attracted to one of the stable

fixed points and locking of two oscillators occurs. Otherwise, a solution with an
n−1
, the
ever-increasing φ takes place. When (ωn − ωn−1 ) is slightly greater then µKn,n−1

dynamics of φ(t) consists of a series of long plateaus (the bottleneck part of the
dynamics of φ) which periodically undergo a rapid slip: a “staircase” solution with
some average growth rate. A period τ over which this phase difference grows by 2π
The appearance of the overbar in δ ω̄n and its lack in δωn±1 is on purpose: while ω̄n represents
the RG solution to the actual running frequency of the crazy oscillator, the frequency correction to
oscillators n − 1 and n + 1 is a renormalization of the parameters that appears in the equations of
motion of these oscillators, but does not yet represent the final value of the running frequency of
those oscillators.

108
can be defined by
Z 2π
τ =

dφ
n−1
(ωn − ωn−1 ) − µKn,n−1
sin φ
2π

= r

(ωn − ωn−1

)3 −

Kn−1
µn,n−1

´2 .

(6.47)

So
φ = θn −θn−1 = (ωn −ωn−1 )

1−

Kn−1
µn,n−1 (ωn − ωn−1 )

¶2
t+oscillating terms. (6.48)

Also notice that that mn−1 θ̇n−1 + mn θ̇n = mn−1 ωn−1 + mn ωn , so
Θ ≡ mn−1 θn−1 + mn θn = (mn−1 ωn−1 + mn ωn )t.

(6.49)

From Eqs. (6.48) and (6.49) we see that mi θ̇i = [mi (ωi + δωi ) + oscillating terms]
where

δωi =

µi,j
(ωj − ωi ) 1 −
mi

1−

Ki,j
µi,j (ωj − ωi )

¶2

,

(6.50)

and i and j are the two oscillators in question (n and n − 1). Note the reciprocity:
mi δωi = −mj δωj .
This was a 2-body system: not surprisingly a fully solvable case. Going back
to the many-body system of a crazy oscillator and the rest of the chain, we will
neglect corrections due to second nearest neighbors and beyond. Then the frequency
corrections are simply the superposition of two 2-body effects. In other words, each
of the nearest neighbors of the crazy oscillator receives a frequency correction given
by Eq. (6.50), and its own frequency correction is given as
δ ω̄n = −

mn+1
mn−1
δωn−1 −
δωn+1 .
mn
mn

(6.51)

If the square roots are expanded to second power in 1/Ω, we recover exactly the same

109
frequency corrections summarized in Eqs. (6.44)-(6.45). Square root formulas of this
type include some of the higher-order effects (such as Ω−ωi in the denominator, rather
then Ω). There may be other terms of this order, for example ωn−1 ωn+1 /Ω2 , which
can not be obtained from a pair of independent 2-oscillator analysis. Therefore, if
used in the full oscillator chain, these expressions form an uncontrolled approximation
at orders higher then O(1/Ω), but empirically we found them to produce a slightly
better match with the simulation data (below) in some cases; we explore this in the
following subsection. Hence, in the numerical renormalization of the chain, we employ
the square root formulas of the type of Eq. (6.50).

6.3.1

Comments on the Crazy Oscillator Step

Now is the appropriate time to reflect on what it means for an oscillator to be crazy.
Notice that θ̇n /ωn ≡ 1 + (K-terms), where the “K-terms” are terms which become
zero when the couplings to the neighbors, Kn−1 and Kn are set to zero. Based on
what we just learned, an oscillator is crazy when these K-terms are ¿ 1. This is just
a definition; we still need to specify what causes an oscillator to possess this property.
Based on this definition, for an oscillator which is crazy, setting the couplings to the
neighbors to zero results in a small correction to θ̇n /ωn , while for an oscillator which is
not crazy, the correction can be substantial. From the reference frame of the oscillator
in question, when the neighbors have significantly different frequencies relative to the
couplings, their pushing and pulling on θn averages out to approximately zero. This
is seen directly from the 2-oscillator analysis, Eqs. (6.46)-(6.50). Hence, setting the
coupling to zero will have a small impact on θ̇n . Thus, it is the frequency difference
relative to the couplings that determines whether an oscillator is crazy. Indeed, the
physics that dictates whether the K-terms are small must clearly be independent of
the reference frame; we can always switch to the frame in which all ω values are large
relative to couplings, so all oscillators appear fast - but not all are crazy! Therefore,
only the frequency differences must enter into the criterion. We define a parameter
rCO as a measure of craziness of an oscillator, or the measure of validity of the

110
corresponding decimation step:

Kn
Kn−1
rCO ≡ Max
µn,n+1 |ωn − ωn+1 | µn,n−1 |ωn − ωn−1 |

(6.52)

Our discussion of the 2-oscillator problem shows that when rCO > 1, the frequencies
of the two oscillators entrain, and for rCO < 1, the corrections due to the interaction
quickly decays, due to the rapid variation of the square root for the small value of
its argument. It is not obvious what the craziness criterion should be in the case of
a crazy oscillator embedded in a many-oscillator case. Numerical experiments, such
as the ones we now present suggest that rCO < 1 is still a reasonable criterion, so we
use it in our implementation of the RG.

Demonstrations
First, consider a system of three oscillators, labelled as oscillators 3, 4, and 5 in
Fig. 6.3 (appears at the end of this chapter). The frequencies of oscillators 3 and 5, as
well as the couplings K3,4 and K4,5 were fixed, while a series of numerical simulations
of Eqs. (6.2) were performed, one for each different value of ω4 ≡ Ω. The resultant
frequencies of each of the oscillators, defined by Eq. (6.40) are measured. This is the
content of Fig. 6.3 (a), along with comparisons with the square-root predictions, such
as Eqs. (6.50),(6.51) and the simpler formula, such as Eq. (6.45).
We see that as Ω becomes large, the average running frequency of oscillator 4, ω̄4
approaches its intrinsic frequency Ω and the average running frequencies of the other
oscillators ω̄i also approach their own intrinsic frequencies ωi . We also see that this
approach is rather rapid and closely follows the square-root predictions.
The detail of the low-Ω region is shown in Fig. 6.3 (b) along with the average
frequencies of all three oscillators, and of the pairs of oscillators 3-4 and 4-5. As Ω is
varied from negative to positive values, first oscillator 3 and 4 form a cluster (point
A), then at even lower Ω all three oscillators entrain to the average frequency (point
B),5 followed by oscillator 3 breaking off from the three-oscillator cluster (point C)

as is expected, since this is the phase-locking situation.

111

Figure 6.3: Top: the set of intrinsic frequencies and couplings. The frequency of
oscillator 4 has been varied. (a) ωi (Ω) for oscillator 3 (green circles), oscillator 4
(red squares) and oscillator 5 (blue diamonds). The square-root approximations,
Eqs. (6.50),(6.51 are shown by the thick curves, the simple formula such as Eq. (6.45)
appears as the dotted curve and the function ω̄ = Ω appears in light grey. (b) Closeup of the low-Ω region. Also shown are the average frequencies of all three oscillators
(thin solid line), of oscillators 3 and 4 (dotted line) and of oscillators 4 and 5 (dashed
line). (c) Same as (a) when oscillators 3, 4, 5 are embedded into a 7-oscillator system
with frequencies and couplings specified in the top plot. (d) ω̄4 (Ω) for three different
functional forms of the interaction. The curve closest to ω̄ = Ω grey line is the sawtooth interaction, the one furthest is the 4-th order polynomial interaction and the
middle curve is the sine interaction, i.e. same as in (a).

112
and finally the 4-5 oscillator cluster breaking up (point D). Moreover, the point of
breakout of oscillator 4 from both of its neighbors (points A and D) in Fig. 6.3 (b) is
close to the point when rCO becomes 1, as the comparison with Fig. 6.3 (a) indicates.
We also embedded the three oscillator system into a larger system with two additional oscillators on each side. The ω̄i (Ω) curves for this situation appears in Fig. 6.3
(c). While the detailed features change somewhat in the low-Ω regions, the rest of the
findings qualitatively remain the same. We note that not all realizations of frequencies and couplings would result in a qualitatively equivalent ω̄i (Ω) plot. For example,
it is possible to come up with such values of parameters for which all three oscillators
are not entrained for any value of Ω. But it appears that it always remains true that
ω̄ of a crazy oscillator approaches Ω for a relatively small values of Ω.
We hypothesized that one of the reasons for this quick approach to the crazy
oscillator regime is the steep decay of the square root correction. The square root
comes from the slowing down as the oscillator passes through the bottleneck of the
washboard potential, see Eqs. (6.46) and (6.47). If this is so, changing the waveform
to increase the duration of this bottleneck should result in a slower approach to the
crazy regime, and changing the waveform to decrease the duration of the bottleneck
should have the opposite effect. To check this hypothesis made numerical simulations
of the 3-oscillator case for three different types of waveforms: sine, the results for
which are already shown, as well as the sawtooth which should decrease the bottleneck slowing down and the periodic quartic polynomial which should enhance the
slowing down. To avoid unnecessary clutter, we define these wave forms graphically
in Fig. 6.4.

The comparison of the results appears in Fig. 6.3 (c). Interestingly, our expectations appear to be confirmed: the sawtooth wave leads to the quickest approach
to the crazy oscillator regime, while the quartic polynomial interaction leads to the
slowest approach. In fact, the analytical form for the period of φ in the case of the

113
Interaction function
1.0

0.5

-6

-4

-2

-0.5

-1.0

Figure 6.4: Triangle wave, sine (middle curve) and quartic polynomial interaction
(outer curve).
sawtooth interaction Saw(φ) is logarithmically diverging:
τSaw

Z 2π
dφ
∆ω 0 1 − κSaw(φ)
1+κ
ln
κ∆ω
1−κ

(6.53)

. For any value of κ this τ is smaller then τSin = ∆ω√2π1−κ2 , the
where κ = µ∆ω

functional form in the case of the sine-wave interaction, so the slowing down of φ and
decrease of its ω̄ relative to Ω is indeed smaller.

6.4

Renormalization Group Implementation

The chain of oscillators is renormalized by successive application of the two decimation steps developed above. These decimation steps are executed numerically on a
list of parameters (mi , ωi , Ki ) representing the chain of oscillators, beginning from
the largest values of K and ω and working down to smaller ones during successive
steps. Initial values of masses in the RG are all set to 1. As already mentioned,
for the purpose of maintaining strong disorder which forms the assumption for the
decimation steps, intrinsic frequencies and couplings in Eqs. (6.2) will be taken from

114
“wide”distributions, i.e., distributions with long tails. The numerical procedure identifies the oscillators with ω and K that lie in a narrow band of magnitudes at the
top of the spectrum of remaining values, and then decimates those oscillators according to the steps defined in Sections 6.2 and 6.3. A band is used, rather than just
selecting the largest values, to improve the efficiency of the code. The width of the
band is chosen to be narrow (1% of the chain) to maintain the descending order in
ω and K for the decimation of nearby oscillators. Oscillators decimated in the crazy
oscillator step form a cluster in the final tally, with the effective mass describing
the size of the cluster and the frequency giving the time-averaged rate at which the
phase of this cluster advances. In decimating these oscillators, we use the squareroot
approximations of the type of Eq. (6.50). As described in Section 6.3, these oscillators are removed from the chain. After the two decimation steps on the band of
oscillators, the remainder of the chain consists of fewer oscillators, which have lower
coupling strengths and frequencies within a narrowed spectrum. The procedure is
then repeated for successively lower bands of K and ω.
During the execution of the RG, we keep track of the values of rSC and rCO . For
most cases in the beginning of the RG procedure, these are not close to 1 because we
work with distributions that exhibit strong disorder. In other words, given a potential
crazy oscillator, in most of the cases the nearby ω and K are small in comparison: this
assumption defines the physical content of the strong disorder. However, it is possible
to have neighbors with parameters which are somewhat lower, yet the difference is not
large enough to guarantee that effects of these neighbors are small. This is especially
important as the spectrum of the remaining K and ω values shrinks. This is where
the importance of rSC and rCO comes in. As we explained in Sections 6.3 and 6.2, we
set set the criteria to rCO < and rSC > 1.
Aside from some rCO and rSC being close to 1, we occsionally encounter another
special case related to violations of the assumptions based on strong disorder. Note
that single oscillators with m = 1 or oscillators representing clusters with m > 1
are normally subjected to the crazy oscillator decimation step at some point in the
RG procedure: most of the bare oscillators either join other oscillators in building

115
clusters, which are eventually decimated via the crazy oscillator step, or they are
decimated as crazy oscillators directly without undergoing a strong-coupling step.
However, there exist rare combinations of frequencies and couplings that will prevent
the decimation from occurring at all. For example, while an oscillator’s frequency
may lie in the executable band, this oscillator may be coupled via a K which is
outside this band (it has a somewhat lower magnitude), but may nevertheless cause
rCO to exceed 1. In general, when the RG reaches the scale of that coupling, it
will perform a strong-coupling decimation step. However, in rare cases this may
not happen due to intervening steps that have modified the neighboring parameters,
rendering that coupling no longer strong according to the rSC criterion. As a result,
the site in question has not been subjected to either decimation step. To address such
loopholes, the algorithm sweeps repeatedly through the entire spectrum of ω and K
values until nothing is left to decimate.

6.4.1

Strain Check

This subsection introduces an additional piece of the RG algorithm which we call
“strain check”. Both of the decimation steps just presented are local: the decision as
to which step to implement and the outcome of each step are based on frequencies
and couplings in the immediate neighborhood of one oscillator or one bond. As
clusters are built up, only two pieces of collective information about the constituents
of the cluster are retained: the average frequency of the constituents’ ω 0 and the total
number of oscillators in that cluster m; all other internal information is forgotten.
However, it is not obvious that only ω 0 and m is enough to accurately predict the
properties of the clusters. We propose an additional measure, the phase strain of a
cluster, not unlike the strain of a material under stress.
It was mentioned in the work of Strogatz and Mirollo [51] that a chain of size
N will phase synchronize, i.e., will allow a solution for θ1 through θN such that

116
θi+1 − θi = const for all i ∈ [1, N ], only when the quantity
Xn =
(ωi − hωi) < Kn ,

(6.54)

i=1

for all 1 ≤ n ≤ N . Here hωi is the average of all bare frequencies of the constituent
oscillators. This conclusion was reached by adding up equations of motion for oscillators 1 through n and making use of the fact that all but the last interaction terms,
Kn sin (θn+1 − θn ) cancel each other out. In order for phase synchronization to happen, the magnitude of the coupling Kn must be greater then or equal to Xn for every
1 ≤ n ≤ N . Notice that Xn is a sum of random numbers, so for larger chains, it is
increasingly likely that Xn will exceed Kn at some oscillator n (see [52] for calculation
of probability that Xn < Kn for all 1 < n < N in the case of all identical Kn = K);
as the size of the cluster increases, the range of values of |Xn | typically grows, and
increasing phase strains |θn+1 − θn | are needed to keep the various parts of the cluster
synchronized. This Xn , which we call “accumulated randomness” serves as a measure of a “strain” of a potential phase-synchronized cluster. When Xn < Kn for all
n, there will be a solution such that all local phase differences lie in [−π/2, +π/2]. If
the strain becomes so large that some |θn+1 − θn | becomes π/2, any further increase
of the strain will cause the cluster to break.
Also note that if two clusters are combined through a strong-coupling decimation
step to form a putative larger cluster, the mean frequency hωi to be used in the
sum in Eq. (6.54) changes, and the condition |Xn | < Kn may become violated for
a bond in the interior of one or other of the component clusters. This means that
the combined set of oscillators are not, in fact, synchronized, and should be broken
into smaller synchronized clusters in the RG procedure. We proceed by making the
assumption that the oscillators form just two synchronized clusters, and identify the
break between the two clusters as where the condition |Xn | < Kn is violated in the
presumed cluster.6 Thus if Xn exceeds Kn at some n, we do not allow the cluster to

Note, it is easy to see that if Xn exceeds Kn only at one place in the cluster, it does not matter
from which end the quantity Xn is computed. If Xn exceeds Kn at more then one place, this only
tells us that the hypothetical cluster would not phase synchronize, but it does not tell where it would

117
form. Instead, we reform the bare oscillators into two different subclusters, such that
the first subcluster is formed out of the bare oscillators 1 through n, and the second
is formed form the remaining bare oscillators. The RG then proceeds as normal. The
new subclusters may be subjected to either of the two decimation steps later in the
sweep, or in the next sweep if more are required.
The strain may appear as incompatible with the idea of the RG because its determination requires the complete information about the constituent oscillators, rather
then some renormalized degree of freedom which hides microscopic details of the original system. However, we propose the strain as a way to check whether additional
information beyond ω 0 and m is necessary. If it turned out that a strain check improves the correspondence of RG and simulation results, it would suggest that we are
on the right track and motivate us to come up with some other method that measures
strain but does not require the knowledge of the oscillators in the original chain. In
our numerical implementation of the RG, the use of the strain check was optional in
order to compare predictions with and without it. We will see in Section 7.2.4 that
the use of the strain check does not seem to provide an obvious improvement in the
correspondence of the RG and simulation data.

break up. In our algorithm, the break was made at the first place where Xn exceeded Kn .

118

Chapter 7
Results
7.1

Testing the validity of the numerical RG

We now switch the discussion from the general RG algorithm applicable to a wide
class of distributions and parameters to specifics of the RG and simulations as they
were used in this work: these include a particular choice of frequency distributions
g(ω) and coupling distributions G(K) as well as other details such as system size, the
set of parameters defining the distributions that are investigated, and cutoff values
in ω and K. Intrinsic frequencies will be assumed to have a symmetric Lorentzian
distribution because they exhibit long tails:
g(ω) =

C1
1 + ω2

(7.1)

with −ωc ≤ ω ≤ ωc . The couplings will be positive and taken from the half Lorentzian
for the same reason:
G(K) =

2µC2
µ + K2

(7.2)

with 0 ≤ K ≤ Kc . The C1 and C2 are normalization constants that become 1/π
when the cutoff values ωc and Kc are infinite (i.e., no cutoff). Notice that due to the
structure of Eqs. (6.2), it is always possible to divide both equations by a constant
such that the width of one of these distributions will be unity (assuming the width of
that distribution is not zero to begin with) - this only changes the timescale of all the
processes in the chain. So, without loss of generality we chose to define g(ω) to have

119
the width of unity and explore the physics by varying the other width parameter µ.

7.1.1

RG Analysis of Lorentzian Chains

The RG procedure is applied to a chain of 106 oscillators. Twelve values of coupling
width µ are studied: µ = 0.025, 0.05, 0.25, 0.5, 0.625, 1.25, 2.5, 3.75, 4.5, 5, 6.25, and
7.5. The set of random numbers defining the particular realizations of the chain were
different in the RG and in the simulations (below) for the statistical comparisons, but
were identical for the real-space comparison described in Section 7.1.3. In order to
facilitate a better comparison with simulations, we chose the same cutoff values for
{ω} and {K} as in the simulations: ω values were chosen to lie in [−100, 100] and K
values were chosen to lie in [0, 100].

7.1.2

Simulations

To better understand the system characteristics and the reliability of the renormalization group, Eqs. (6.1) are also integrated numerically. The numerical method is
a variable stepping Runge-Kutta algorithm. Systems of N = 10000 oscillators are
solved with N intrinsic frequencies chosen randomly from the symmetric Lorentzian
given by Eq. (7.1). Similarly, N − 1 coupling constants are randomly selected from
the half Lorentzian given by Eq. (7.2). For both distributions the same cutoff value
is used: ωc = 100 and Kc = 100; increasing cutoff values increases the time needed
for simulations. The same twelve values of coupling width µ that are used in the RG
are studied in the simulations: µ = 0.025, 0.05, 0.25, 0.5, 0.625, 1.25, 2.5, 3.75, 4.5,
5, 6.25, and 7.5. Results for each value of µ are averaged over 100 realizations of the
intrinsic frequencies and coupling constants. To help facilitate comparisons between
different distribution widths the same 100 random number seeds are used for each
coupling width.
The simulations are integrated for a relatively long time of t = 10000. Running
phases at each site are recorded at regular intervals and used to calculate average

120
frequencies over a time T according to
ω̄i (T ) =

θi (t0 + T ) − θi (t0 )

(7.3)

with T = t − t0 [compare with the theoretical definition of ω̄ given by Eq. (6.40)].
The time t0 is chosen to eliminate transients: increasing values are tried until average frequencies remain unchanged. In these simulations, the value of 2π/T sets
the resolution limit in ω̄i to distinguish two neighboring frequency clusters. A group
of oscillators is determined to be members of a synchronized cluster only if all the
members (i = n, n + 1, . . . , n + m + 1) have the same value of ω̄i as the neighboring
sites within some tolerance |ω̄` − ω̄i | ≤ η (` = i − 1 and i + 1). Since T −1 is O (10−4 )
the tolerance η is set to 10−3 .

7.1.3

Comparison of RG and Numerics

Insight into the validity of this RG can be found by careful comparison of the frequency cluster structure that it predicts with results from numerical integration of
Eqs. (6.2). If the renormalization method is wrong there is expected to be little
agreement between the two. Figure 7.1 plots these predictions over a representative
sample of 150 oscillators from chains of N = 10000 at 3 different coupling distribution
widths. For all 3 values of µ, the RG predicts clusters with large frequency oscillators
interspersed. As expected, the characteristic size of these clusters grows with µ. Also
shown are time-averaged frequencies, Eq. (7.3), calculated from numerical solutions
of Eqs. (6.2). There is excellent agreement of the cluster sizes, locations, and frequencies between the two approaches. This provides some confidence that the essential
mechanisms are captured by the RG procedure.
Analysis of these and similar plots, shows that the agreement of the cluster frequencies becomes more accurate with increasing cluster size. We do in fact expect the
dynamics of most of the oscillators in large cluster to be very similar to the dynamics
if that cluster was isolated from the rest of the chain. The reason for this is the
smaller boundary-to-bulk ratio of larger clusters, which is more apparent if Eq. (6.2)

difference

simulated
dynamic
dynamic
frequencies
frequencies

121
4700 4800 4900 5000 5100

4700 4800 4900 5000 5100

(a)

4700 4800 4900 5000 5100

(b)

(d)

(c)

(e)

(f)

50
25
-25
-50
-1
-3
-5

(g)

4890 4900 4910 4920

(h)

4890 4900 4910 4920

(i)

0.5
-0.5

4890 4900 4910 4920

Lattice Site

Figure 7.1: Cluster structure of a small section of a 10000-oscillator chain. The three
columns correspond to µ = 1.25, µ = 3.75 and µ = 7.5 respectively. The top row
displays 500 time-averaged running frequencies ω̄ from simulations of Eq. (6.2). The
middle row zooms in on a smaller section of 50 oscillators and compares ω̄ from
simulations (solid squares) and from the RG (open circles). The bottom row displays
the diﬀerence between the frequencies ω̄ of the RG and the simulations for the same
50 oscillators as in the middle row.
is rewritten as
θ̇i = ωi +

Ki−1
Ki
sin (θi−1 − θi ) +
sin (θi+1 − θi ).
mi
mi

(7.4)

As clusters become larger, the importance of the coupling to the rest of the chain is
weighted down by a 1/m factor. In the crazy oscillator decimation, the boundaryto-bulk ratio is manifested by the appearance of the 1/m factor in the frequency
correction. Notice that in time-averaging Eq. (6.2), the inﬂuence of 1/m suggests
that
ω̄i ≡ θ̇i t ≈ ωi ,

(7.5)

but ωi is given by the mass-weighted average of the constituent bare ω values. Also,
when a cluster is built up by a strong-coupling decimation step, any frequency renor-

122
malizations of the component oscillators by prior fast oscillation decimation steps
cancel when the weighted sum of the frequencies is calculated to give the frequency
of the new cluster. Thus the cluster frequencies ω̄i are given approximately by the
average of the constituent bare ω values, plus the small correction due to the crazy
oscillator decimation step that yields the final cluster. This is what we see in comparison of RG with simulations: the ω̄ of clusters in simulations is accurately predicted
by the RG predictions, which are mostly given as by the mass-weighted average of
the constituent bare ω values (plus the small correction due to the crazy oscillator
decimation step).
The agreement with exact numerics in a small chain gave us reassurance that the
RG method works, at least in some parameter regimes. Investigating the frequency
distribution of clusters of a given size provides further insight into the validity of this
RG procedure, but this discussion is postponed until Section 7.2.4.

7.2

Physics of the Unsynchronized Phase

We are now ready to turn to the prime goal of this work - understanding the behavior
of the random oscillator chain, and its universal aspects. We first use the results
of the RG and the simulations to discuss the statistical properties of the frequency
clustering in the one-dimensional oscillator chain with strong randomness.

7.2.1

Cluster Size Distribution

The central quantity for understanding the behavior of the random oscillator chain is
the distribution of cluster sizes. In Fig. 7.2 we show the distribution of cluster sizes
for several widths µ of the coupling distribution from the RG and simulations. The
RG was performed on one single realization of 106 oscillators, while the simulations
were performed on 100 realizations of 104 oscillator chains. Since the characteristic
cluster size is much smaller then 104 , we combine results from all 100 realizations and
analyze them as a single chain of size 106 . As can be seen in Fig. 7.2, good agreement
exists between the RG and the simulations. The apparent discrepancies at the largest

123
cluster sizes are due to statistical fluctuations resulting from the small number of such
clusters present in the 106 oscillator ensemble.

Relative number

µ=7.5

µ=5
−2

µ=3.75

µ=2.5
µ=1.25

−4

Cluster size

Figure 7.2: Number of clusters of a given size versus the size. Comparison between
simulations of Eq. (6.2) (solid squares), the RG without strain check (open circles)
and the RG with the strain check (diamonds). The curves are spread out for clarity
by multiplying each successive data set as µ decreases by 0.1.
For a given value of the coupling width parameter µ, the probability of finding a
cluster is predicted to fall off rapidly with increasing cluster size. The linear scaling
in this semilog plot suggests that distribution of cluster sizes has the form of an
exponential: P (n) ∝ exp (−n/ξ), where P (n) is the number of clusters of a given
size, for example. The characteristic length ξ was extracted by making a linear fit
to the data (not shown in Fig. 7.2) between twice the estimated characteristic size ξ
(so it is an iterative process), and the cluster size at which there are 20 occurrences

124
of clusters of that given size.1 The following section focuses on the dependence of ξ
with µ.
We also investigated the effect of the larger system size. In Fig. 7.3 we show
comparison of the RG for a system of a one million and ten million oscillators. Aside
from continuing the exponential trend for larger cluster sizes, there does not seem to
be any significant extra information conveyed by the larger system size.
−1

−1

P(n) / (total number of oscillators)

(b)

(a)
−2

−2

−3

−4

−5

−6

−7

10
10

Cluster size

Figure 7.3: Number of clusters of a given size P (n), normalized by the system size
for µ = 7.5 at two system sizes: 106 oscillators (open symbols) and 107 oscillators
(smaller filled symbols). (a) RG without strack check and (b) RG with the strain
check.
The values of ξ were obtained as described in the previous paragraph. We find that
in the system with 106 oscillators, ξ = 3.28 and ξ = 3.05 without and with the strain
check respectively, while in the system with 107 oscillators, ξ = 3.26 and ξ = 2.97
without and with the strain check respectively.
The distributions P (n) are not strictly exponential, but they appear to be asymp1

since the standard deviation is expected to be of order of the square root of the number of
occurrences.

125
totically so, as seen from the plot of P (n) for a system of 107 oscillators and the fit
which was made between twice the characteristic length (7) and the cluster size 35 at
which there are 27 clusters of that size: see the low-n part in Fig. 7.4 (a). Also shown
in Fig. 7.4 (b) is the plot of slopes of partial fits, i.e., fits to P (n) between cluster size
m and cluster size 35. Aside from the low-statistics scatter at the end, we see that
the variation in this data is much smaller then the average, and the data appears to
be constant, indicating asymptotically exponential P (n).
−0.1

(b)

(a)

Slopes of partial fits

P(n)

−0.2

−0.3

Cluster size

−0.4

Figure 7.4: (a) P(n) along with the fit between cluster size 7 and 35 as described in
the text. (b) The slopes of partial fits between cluster size m and cluster size 35.

7.2.2

Dependence of ξ on Disorder Strength

The width of the coupling distribution, µ is the only control parameter that characterizes the chain. Also, as we showed in the previous section, distributions of cluster sizes
are asymptotically exponential, and hence they are characterized by a single parameter ξ. Hence, the function ξ(µ) is a comprehensive quantity that characterizes the
physics of the chain in the unsynchronized regime. We plot this function in Fig. 7.5.
We find good agreement between simulation (solid dots) and RG (open circles) until
µ ≈ 2.5, after which the two results begin to deviate from each other. The small,

126

3.5

(a)
2.5

1.5

0.5

(b)

0.5

−2

−1

Figure 7.5: Characteristic cluster length ξ given by the inverse of the fits to data in
Fig. 7.2. (a): Results of simulation (filled squares), the numerical RG of Section 6.1
(open circles) and an enhanced numerical RG with strain check (diamonds) described
in Section 7.2.2. The solid line is a simple prediction ξ(µ) = −1/ log(1 − ℘(µ)) where
℘(µ) comes from Eq. (7.8). The dashed line is a linear fit to the simulation data
(solid squares) for µ > 0.625 and has a slope of 0.461; similar fits to the circles give
the slope 0.477 and to the diamonds the slope of 0.478. (b): Log-log plot of the same
data excluding the analytic prediction.

127
0.65

0.6

0.55

Running slope

0.5

0.45

0.4

0.35

0.3

0.25

0.5

1.5

2.5

3.5

4.5

5.5

6.5

7.5

ln(µlower)

Figure 7.6: Slopes of partial fits of ξ versus ln µ between ln µlower and ln 7.5. Solid
dots: simulation; open circles: RG without strain check; diamonds: RG with the
strain check.
but growing discrepancy beyond that µ suggests two possibilities: (i) existence of an
additional mechanism in the physics of cluster formation that becomes important at
larger µ and is not taken into account by the decimation steps of the RG procedure;
(ii) corrections are needed to the steps already included in the RG procedure. The
ξ(µ) obtained with the version of RG which included the strain check is shown in
diamonds in Fig. 7.5. All three sets of data, when plotted on a log-log scale follow
approximately straight lines for µ > 0.5, so ξ(µ) appears to become a power law in
this range. The straight-line fits to the data were made for µ between 7.5 and 0.625.
The error bars in the values of these exponents can be estimated from the scatter in
the partial fits in Fig. 7.6 (we chose the range 0.625 < µ < 7.5 for doing this). This
leads to the following estimate of ξ: 0.47 ± 0.01 for simulations, 0.48 ± 0.02 for RG
with and without the strain check (when rounded to two digits).

128

7.2.3

Approximate Analytical Argument

A simple estimate of characteristic cluster sizes that reproduces the exponential form
of P (n) and also puts an upper bound on the characteristic cluster size can be made by
noting the role of weak couplings in limiting cluster sizes.2 Indeed, that small clusters
are likely to form in between weak couplings (also between the crazy oscillators,
but these are padded by a weak coupling on each side, so are automatically taken
into account), and since the distribution of separations between such weak couplings
follow Poisson statistics the exponential distribution should be expected. Using this
hypothesis the probability of finding a cluster of size n (i.e., a cluster delimited by
two weak couplings n units apart) is P (n, ℘) = C(1 − ℘)n = Ce−n/ξ where ξ =
−1/ ln (1 − ℘) and ℘ is a probability that any randomly chosen coupling is weak. We
make an estimate of this ℘ based on bare distributions of ω and K. Similarly to the
RG, we define a coupling to be weak if 2K < |ωl − ωr |. For each such bond, the
probability that the required ratio is > 1 is given by

Z ∞ µZ ∞

G(K) dK

g(ω)g(ω + δ) dω ×

℘=2

ÃZ
δ/2

dδ.

(7.6)

−∞

The integral of g(ω)g(ω + δ) is the probability that for each of the two oscillators
connected by this bond, ωl − ωr = δ. For each such frequency difference, all the
bonds with K < |δ|/2 are considered weak. Finally, all possible values of δ are
integrated in the outer integral. After performing the inner integrals over ω and K
we obtain
℘(µ) = 3

Z ∞

2π
arctan
4 + δ2

2µ

dδ.

(7.7)

To handle this integral we differentiate ℘ with respect to µ under the integral sign,
eliminating arctan and permitting the resulting expression to be integrated in δ in

The exponential dependence of the correlation function, which is closely related to distribution
of cluster sizes, is actually rather familiar from equilibrium statistical mechanics. Examples include
the Ising model and XY model that display exponentially decaying correlations [102]. This can be
shown using standard techniques, such as transfer matrices, or by performing a high-temperature
expansion [50]. However, underpinning these explanations of the exponentially decaying correlations
is a partition function, which is inapplicable in our non-Hamiltonian model, so we are forced to seek
a different explanation

129
closed form. Finally, by reintegrating with respect to µ and redefining variables once
more we find

Z ∞
℘(µ) =

4 ln (x2 /4)
dx.
2µ π (x − 4)

(7.8)

This result is compared with both the RG predictions and numerical solutions in
Fig. 7.5. The expression gives a good description of the data for µ → 0, but overestimates the characteristic cluster size for larger values of µ. Of course such simple
estimate for ξ is but an upper bound; it does not take into account any information
about the constituent ω or K. The renormalization group method, which does not
treat all scales of frequencies and couplings at once, clearly provides better predictions then an estimation based on bare distributions. The noteworthy feature of the
RG procedure is that the process of building up clusters by successive applications
of very simple decimation steps leads to good predictions of structures which are far
from trivial!

7.2.4

Cluster Frequency Distributions

Investigating the frequency distribution of clusters of a given size provides further insight into the physics and the validity of this RG procedure. Distributions in Fig. 7.7,
serve as representative examples.3

The RG predictions are compared with the re-

sults from simulations of Eqs. (6.1) for several clusters sizes and parameters µ. Notice
the excellent agreement for the larger values of m and also for all m at small values
of µ. The comparisons in Fig. 7.7 do show some underestimation of the number of
smallest sized clusters (m = 1 and m = 2) near zero frequency. The region of the
greatest discrepancy, small m and ω, is in fact where the crazy oscillator decimation
is expected to be least accurate, since by the time the RG procedure reaches small
frequencies, the distributions of remaining frequencies and couplings are no longer
expected to be wide, so the approximations based on strong randomness are not as
accurate. However the disagreement decreases rapidly with increasing m, because as
A reader is invited to see [101] for a different, but complimentary way of plotting these distributions.

130

µ = 1.25

m=1

1000

500

2000

500

200

−10

300
−5

2400

1500

1200

700

1500

−5

700
300

−5
1400

−5

700
700

700

−5

m=4

−10

1000

2000

m=3

µ = 7.5

4000

−10
4000
m=2

µ = 3.75

−2

1000

500

300
−2

−4

−2

700

m=5

−2

200
−4

Figure 7.7: Frequency distributions. The y-axis displays the number of clusters of
a given size m in a frequency bin of width ∆ω = 0.1 at frequency specified on the
x-axis.We display distributions for m = 1, 2, 3, 4 at µ = 1.25, 3.75 and 7.5. Thick
curve: numerical simulation; thin red curve: numerical RG without the use of the
strain check.

131

µ = 1.25

m=1

1000

500

2000

500

200

2000

m=3

−5

−10

300

1200

700

1500

−5

700
300

−5
1400

700

−2

−5

300
−2

−4

1000

700

500

300

m=5

700

−10

700

1500

700

1400

2400

−5

m=4

µ = 7.5

4000

−10
4000
m=2

µ = 3.75

−2

−4

−2

Figure 7.8: Frequency distributions. The y-axis displays the number of clusters of
a given size m in a frequency bin of width ∆ω = 0.1 at frequency specified on the
x-axis.We display distributions for m = 1, 2, 3, 4 at µ = 1.25, 3.75 and 7.5. Thick
curve: numerical simulation; thin green curve: numerical RG with the use of the
strain check.

132
described above, the frequencies of larger clusters are insensitive to the corrections
from neighbors: for large m the running frequency is mostly given by the average of
intrinsic values of ω. We also note in passing that the disagreement decreases with
decreasing µ, but refrain from proposing an explanation for this.
Another remark is regarding the shape of distributions for m = 1 oscillators: they
do not fit pure Lorentizian at all values of ω, however at very large ω the tails of
RG and simulation results match the tails of the bare frequency distribution. This is
expected because these tails are comprised of very high frequency oscillators which are
typically decimated out as crazy oscillators to form one-oscillator “clusters.” Because
of their high bare ω, the frequency correction from the neighbors is typically negligible,
so their renormalized frequencies equal bare frequencies. An interesting observation
however (not shown in Fig. 7.7), is that the frequency interval over which there is a
good match between RG and simulation for these m = 1 oscillators extends to lower
frequencies then where the RG and simulation data match with the bare Lorentzian,
providing yet another indication for the validity of the frequency adjustments included
in the fast oscillator decimation steps.
The use of strain check in the RG appears to increase the number of low-frequency
clusters for m not too large and for all sampled µ, as can be seen from Fig. 7.8.
Although not obvious from these figures, closer inspection reveals that the tails are
not significantly affected - the increase is mostly at low frequencies. This net increase
in the number of small-m clusters is compensated at larger m - this is not shown in
Fig. 7.8 but can be seen from closer inspection of large-m distributions and can also
be seen from statistics of cluster sizes, Fig. 7.2 (notice the diamonds slightly above the
circles for small m and below the circles for larger m). Indeed, the strain check acts
to mostly break up large clusters, therefore decreasing their number and increasing
the number of smaller ones, but it is not trivial to see why this increase affects mostly
low frequencies.
We leave this section with a reminder that a more refined measure of a cluster
strain or other improvement in the RG strategy, which would reduce differences
between all RG results and simulation results, is due to be invented in the future.

133

Chapter 8
Concluding Remarks
In this work we investigated the collective behavior of a one-dimensional chain of
coupled nonlinear oscillators with random frequencies and nearest-neighbor couplings.
For this purpose we developed a real space renormalization group approach, which is
expected to be reliable in analyzing chains with large disorder. As we report above,
our RG approach did an excellent job in capturing both the dynamics of individual
oscillators, and the large scale behavior of the chain.
The disordered oscillator chain has the potential of establishing macroscopic order,
which is exhibited by a finite fraction of all oscillators in the system moving in accord
with the same frequency. Indeed in systems with higher connectivities, such global
synchronization may occur, but in the limit of short range interactions, the macroscopic order is stymied by fluctuations. This dynamical behavior is reminiscent of
critical phenomena in equilibrium statistical mechanics, where for each system there
exists a lower critical dimension, and only dimensionality higher than that allows
macroscopic ordering. In our analysis, we concentrated on the quantities that most
directly capture the collective aspects of the chain’s behavior, and therefore also the
competition between the interactions which seek to establish a synchronized phase,
and the disorder-induced fluctuations which destroy it. This physics is fully contained
in the cluster-size distribution, P (n) ∼ e−n/ξ(µ) , which we find as a function of the
interaction-strength tuning parameter, µ.
The behavior of the average cluster size, ξ(µ), acts as an effective correlation
length for the oscillator chain, where we once more think of our system in terms of

134
Unsynchronized
fixed point

Fully synchronized
fixed point

Figure 8.1: Fixed points and flow diagram of the oscillator chain. The short-range
oscillator chain has only two fixed points: the unstable fully synchronized fixed point
at infinite interactions (or zero disorder) µ → ∞, and the stable unsynchronized fixed
point. The cross over flow between the two points is associated with a correlation
length, captured by the average cluster size, ξ(µ).
the theory of critical phenomena. Since our oscillator chain is below its lower critical
dimension, its phase diagram resembles that of the one-dimensional Ising model, with
µ serving the role of inverse temperature 1/T . The chain has an ordered phase, where
all oscillators are synchronized, which is reached in the limit µ → ∞ (analogous to
the Ferromagnetic phase at T → 0). This phase, however, is described by an unstable
fixed point. If an RG flow can be extended to the entire parameter space, the fullsynchronization fixed point flows under RG to the stable fixed point found at the
noninteracting point, µ = 0, which describes the unsynchronized phase (analogous to
the paramagnetic fixed point at T → ∞). The putative flow diagram of the chain is
shown in Fig. 8.1. The flow from the full-synchronization point to the unsynchronized
µ = 0 point is indeed associated with a length scale, whose meaning is the coarsegraining length at which the chain seems fully unsynchronized, i.e., consisting of single
free oscillators. This length is the size of synchronized clusters in the original chain,
ξ(µ). Since ξ(µ) serves as an effective correlation length, it may diverge in a universal
fashion in the limit µ → ∞. While our RG method can only access moderate values
of µ, we cautiously argue that we can access the critical region, as can be seen in
Fig. 7.2. From our results we obtain the critical exponent ν of the synchronization
transition:
ξ(µ) ∼ µν

ν ≃ 0.48 ± 0.02

(8.1)

Note that this behavior is not analogous to the divergence of the Ising chain correlation
length ξ ∼ ln T1 .
The method by which we obtained the above results constitutes the main achieve-

135
ment of this work. We developed a real space renormalization group scheme that
successfully predicts the detailed behavior of a nonlinear oscillator chain with shortrange interactions. We have implemented the RG scheme numerically on a chain of
106 oscillators and compared its results with exact simulations of the model equations,
Eqs. (6.1), in chapter 7. The agreement is very good, and for some characteristics,
such as the cluster frequency distributions at low frequencies, it improves with larger
cluster sizes. This indicates particularly that our local RG scheme successfully captures the nontrivial collective features of the chain expressed in large synchronized
clusters.
One practical implication of our analysis may be found when considering coupled
laser arrays, where only limited connectivity can be achieved. Diode lasers, for example, are manufactured in one-dimensional bars that allow for evanescent coupling
through overlapping electric fields. Due to the exponential falloff in the electric field
envelope with distance, the coupling between well-behaved diodes can only be short
ranged. The tools developed here could be used to study the extent of coherence that
such systems exhibit. A potentially promising area of current research is constructing grids of diodes coupled in 2-dimensions by stacking 1-dimensional bars. Our RG
technique could possibly also be extended to such a system, and help clarify both its
properties, and whether it can establish macroscopic coherence.
Our work suggests several other directions for future research. First, we may
ask how can one improve the applicability of our RG scheme to larger values of the
interaction parameter µ. This may be pursued by focusing on the concept of intracluster strain, which is roughly the variance of phase differences between neighboring
oscillators within a cluster. As interaction increases, cluster sizes grow, and bigger
strain is needed to accommodate the spread of frequencies within a cluster (cf. Sec.
7.2.2). In our current scheme we made the first attempt to take this into account with
the strain check based on the argument of [51], a condition expressed by Eq. (6.54).
This led to some improvement - for example, ξ(µ) was modified in the right direction,
but this modification was overestimated. The strain check may be augmented in the
future.

136
A second natural question is the degree of universality in the oscillator chain. Our
result for the correlation length critical exponent ν, Eq. (8.1), was obtained using
Lorentzian distributions. It is possible that the exponent depends on the nature of
the tails of the distributions used; Lorentzians, for instance, do not have a variance.
We intend to investigate the universality of ν alongside investigating the applicability
of our RG scheme to distributions with narrower tails. This work is already underway
[103].
A more general question concerns the development of an RG scheme and the analysis of synchronization in short-range networks at higher dimensions, and especially
at two dimensions, which is very close to the lower critical dimension: Refs. [104, 105]
suggest that in two dimensions there may be a transition to a macroscopically synchronized phase, or that d = 2 is the lower critical dimension giving rise to interesting
anomalous behavior (see also [45]). The development of a decimation procedure at
dimensions higher than one is quite challenging because of the higher connectivity
of the system, and the possibility that real-space local decimation steps change the
geometry of the system.

137

Appendix A
Relating Beam Theory to Duffing
Equation
In this appendix, the parameters of the Duffing equation (except damping) are extracted from the underlying equation of motion of a doublyclamped beam. This is
useful to obtain an estimate of the realistic values of the parameters that go into the
Duffing equation and also to see explicit dependence of these parameters in the two
regimes: the tension-dominated regime and the stiffness-dominated regime. In the
process of preparation of this appendix, I learned that the review article [5] contains a
similar discussion. I chose to include this appendix in any case, because it was part of
the work I have done while collaborating with experimentalists. Also, figures are included here which do not appear in [5]. A different method (Galerkin approximation)
has recently been used for similar calculations in [4].
The beam will be described by a 1-dimensional displacement function Y (x, t). We
start with the following equation describing the dynamics of the beam:
Z µ ¶2 # 2
Y A l ∂y
∂ y
∂ 2y
∂4y
dx
−ρA
Y I 4 − τ0 1 +
∂x
2l 0 ∂x
∂x2
∂t2

(A.1)

The integral term corresponds to increasing tension when the beam is deformed, i.e.,

138
the tension depending on the shape of the beam. Let us first nondimentionalize:
y = ry 0

(A.2)

x = l 0 x0

(A.3)

t = t0

ρAl02
τ0

(A.4)

The resulting dimensionless equation (dropping 0 ) is:
Z l/l0 µ ¶2 # 2
∂y
∂ y
∂ 2y
∂ 4y
ξ 4 − 1+ν
dx
∂x
∂x
∂x2
∂t2

(A.5)

where
YI
l02 τ0
Y A r2
ν =
2τ0 ll0
ξ =

(A.6)
(A.7)

For a double-clamped beam, l remains strictly equal to l0 (and to a good accuracy for
a cantilever as well, although technically it can stretch). Also, for a doublyclamped
beam the boundary conditions are
∂y ¯¯
∂y ¯¯
y(0, t) = y(1, t) =
= 0.
∂x ¯(0,t)
∂x ¯(1,t)

(A.8)

For ν = 0 the equation is linear with the general solution
y0 (x, t) =

Yn (x)e−iωn t ,

(A.9)

n=1

where Yn are eigenfunctions of the operator L̂:
L̂Yn (x) = ωn2 Yn (x),
d4
d2
L̂ = ξ 4 − 2 ,
dx
dx

(A.10)
(A.11)

139
with boundary conditions
Y (0) = Y (1) = Y 0 (0) = Y 0 (1) = 0.

(A.12)

It is straightforward to show that with these boundary conditions, L̂ is Hermitian,
i.e.,

Z 1

Z 1
Yn (x)L̂Ym (x),.x =

Ym (x)L̂Yn (x),.x.

(A.13)

Therefore, its eigenmodes form a complete set, so a nonlinear solution can be expanded in terms of them
y(x, t) =

a(t)Yn (x).

(A.14)

n=1

Analogously to the time-dependent perturbation theory in Quantum Mechanics the
nonlinearity will spread solution among the linear modes, even if initially the initial
condition is fully in one single mode. Nevertheless, we will make an approximation
based on the physical grounds that when the amplitude is sufficiently small, the
system is approximately in the ground state, i.e., y(x, t) ≈ a(t)Y0,1 (x), and we would
like to know the effect of nonlinearity on the timedependence of y, i.e., on a1 . In this
spirit, we substitute a(t)Y0,1 (x) into Eq. (A.5), multiply by Y0,1 (x) and integrate out
the x-dependence. The resulting equation for a(t) is
ä = −ω 2 a(t) − γa3 (t),
where

(A.15)

R 1 ³ ∂ 4 Y1 ∂ 2 Y1 ´
ξ ∂x4 − ∂x2 Y1 (x) dx
ω2 =
R1 2
dx
0 1

and

³R
γ = −µ

1 ∂ 2 Y1
Y (x) dx
0 ∂x2 1

´ ³R ¡ ¢
1 ∂Y 2

Y12 dx

∂x

(A.16)

(A.17)

Actually, Eq. (A.16) is identical to Eq. (A.10) with n = 1, so in the approximation
used, the coefficient ω does not depend on the presence of nonlinearity. Eq. (A.16)
can also be seen as a variational solution to ω. However, we do not need to employ

140
the variational method (which by its nature gives an approximate values) because it
is best to solve the eigenproblem defined by Eqs. (A.10)-(A.12) directly. As an ODE,
Eq. (A.10) has the following general solution
Yn (x) = an eλn x + bn e−λn x + cn cos (λ′n x) + dn sin (λ′n x),

(A.18)

where

λn =

1 + (1 + 4ξωn2 )1/2
2ξ

(A.19)

λ′n =

−1 + (1 + 4ξωn2 )1/2
2ξ

(A.20)

However, Eq. (A.10) is a boundary-value problem which needs to satisfy the boundary
conditions given by Eq. (A.12). These conditions will determine the eignevalues
{λn , λ′n }, as well as the coefficients {an , bn , cn , dn }. The doublyclamped boundary
conditions lead to the following equation (dropping the subscripts n for a moment):


   
 λ
 e
e−λ
cos λ′
cosλ′   b   0 
  =  
   
 c   0 
 λ
−λ
   
−λ
λe −λe
−λ sin λ2 λ2 cos λ2

(A.21)

For a nontrivial solution we demand that the matrix on the left hand side is noninvertible by setting its determinant to 0. This gives an implicit equation for the set
{ωn (ξ)}:
f (ω, ξ) = 4

ωn2
[1 − cos λ′n cosh λn ] + sin λ′n sinh λn = 0,

(A.22)

where λ′n and λn are given by Eqs. (A.19)-(A.20). For a given value of ξ, there are
an infinite number of zeros of the function f at a discrete set of ωn . Using this, we
can substitute the resulting Y1 (x) into Eq. (A.17) to obtain γ. See Fig. A.1.

141

Figure A.1: Top: ω12 (ξ)/ξ. Notice the crossover from the the string-dominated regime
at lower ξ to the stiffness-dominated regime at larger ξ. Bottom, −γ/µ.

142

Appendix B
Properties of Noise in the
Amplitude Equations
Referring back to definitions in 2.1.1, we note that
31/2 ³ η ´3/2 Ω τ +2π/Ω ˜
NQ (T ) =
f (ξ) sin (Ωξ) dξ
π τ
31/2 ³ η ´3/2 Ω 2ηT /²+2π/Ω ˜
f (ξ) sin (Ωξ) dξ.
π 2ηT /²

(B.1)
(B.2)

Similarly
31/2 ³ η ´3/2 Ω
NP (T ) =

Z 2ηT /²+2π/Ω

f˜(ξ) cos (Ωξ) dξ.

(B.3)

2ηT /²

We compute correlation functions,

hNQ (T1 )NQ (T2 )i =
µ ¶2 Z 2ηT1 /²+2π/Ω
3 ³ η ´3 Ω

2ηT1 /²

dξ

Z 2ηT2 /²+2π/Ω

dξ 00 hf˜(ξ 0 )f˜(ξ 00 )i sin (Ωξ 0 ) sin (Ωξ 00 ).

2ηT2 /²

(B.4)

143
Using the correlation property of f˜, Eq. (2.10) we have

hNQ (T1 )NQ (T2 )i =
µ ¶2 Z 2ηT1 /²+2π/Ω
Z 2ηT2 /²+2π/Ω
3 0 ³ η ´3 Ω
dξ
dξ 00 δ(ξ 0 − ξ 00 ) sin (Ωξ 0 ) sin (Ωξ 00 ).
2ηT1 /²
2ηT2 /²
(B.5)
We make a simple change of variables y = Ωξ and also the property of Delta functions
δ(ay) = a1 δ(y). Then
hNQ (T1 )NQ (T2 )i =

3 B 0 ³ η ´3
2Ω ²

µ ¶2 Z L1 +2π
Z L2 +2π
dy 0
dy 00 δ(y 0 − y 00 ) sin (y 0 ) sin (y 00 ),
L1
L2
(B.6)

where Li = 2ηTi Ω/². Similarly,
µ ¶2 Z L1 +2π
Z L2 +2π
dy
dy 00 δ(y 0 − y 00 ), cos (y 0 ) cos (y 00 )
L1
L2
(B.7)
µ ¶2 Z L1 +2π
Z L2 +2π
3B η
dy 0
dy 00 δ(y 0 − y 00 ). cos (y 0 ) sin (y 00 )
hNP (T1 )NQ (T2 )i =
2Ω ²
L1
L2
(B.8)

3 B 0 ³ η ´3
hNP (T1 )NP (T2 )i =
2Ω ²

Let us investigate the properties of
Z L1 +2π
Is =

dy
L1

Z L2 +2π

dy 00 δ(y 0 − y 00 ) sin (y 0 ) sin (y 00 ).

(B.9)

The value of this integral is zero unless the square of integration includes the y 0 = y 00
diagonal line, in which case it is equal to
Is = π −
|∆| + sin (2L1 ) − sin (2(L1 + |∆|)) ,

(B.10)

144
where ∆ = L2 − L1 . In other words,
µ ¶2
3 B 0 ³ η ´3 Ω
hNQ (T1 )NQ (T2 )i =
2Ω ²
 π − 1 £ 2ηΩ ∆T + 1 sin ¡ 4ηΩ T ¢ − 1 sin ¡4 ¡ ηΩ T + ηΩ ∆T ¢¢¤ |∆| ≤ π² ,
² 1
ηΩ
π²
|∆| > ηΩ
(B.11)
where ∆T = T2 − T1 . This finding is very peculiar. First, the correlation function is
not a delta function in ∆T , so the noise terms N are not strictly white, and second it
is dependent not only on the time difference, but also on the time itself! The integral
π²
π²
of this autocorrelation function over − ηΩ
< ∆T < ηΩ
is

3 B 0 ³ η ´3
2Ω ²

µ ¶2 # · 2
¶¸
π ²
4ηΩT1
π²
sin
ηΩ
2ηΩ

(B.12)

However, the support of the autocorrelation function of NQ is very narrow (compared
to the height) and the fast time dependence of the area averages out to zero on the
long timescale, so we will make an approximation
µ ¶2 #
π2²
3 B 0 ³ η ´3 Ω
δ(∆T )
hNQ (T1 )NQ (T2 )i ≈
2Ω ²
ηΩ
3 γ ³ η ´2
Bδ(∆T )
2 ω05 ²
3γη 2
B δ(∆T ),
8ω03 Γ2

(B.13)

with the identical result for hNP (T1 )NP (T2 )i. As for the cross-correlation functions,
π²
π²
it is odd in ∆T and so its integral over − ηΩ
< ∆T < ηΩ
zero. Thus we consider

hNQ (T1 )NP (T2 )i = hNP (T1 )NQ (T2 )i = 0.

145

Appendix C
Features in the Pattern of Classical
Trajectories
The goal here is simply to provide examples of some of the interesting features observed in the bistable regime which were not mentioned in the main body of the
thesis.

C.1

Shadow regions

Shadows can be seen in the pattern of escape trajectories terminated at the separatrix.
Because there is nothing fundamentally special about the separatrix, there is nothing
special about shadows; they can of course be reached from the attracting fixed point
by following a trajectory which leaves the basin and then re-enters it and eventually
wanders into the shadow region. They are important however, when one is concerned
about the distribution of exit locations along the separatrix: probability of exiting
the basin through a given point on the separatrix, which has been discussed in Maier
and Stein’s long SIAM paper [74].
Here, a specific example will be considered. Using the notation of my thesis, we
will choose η = 0.18 and F ≈ 0.164 which corresponds to the 10% of the hysteresis and
will look at transitions from the lower amplitude branch (the down → up transitions).
Consider first a set of 50 trajectories that escape from the attracting fixed point
depicted in Fig. C.1. All these trajectories satisfy the Hamilton’s equation with energy
0, but they differ in trajectory parameters (see Section 3.1 of the Thesis ). The range

146
of trajectory parameters in Fig. C.1 is [0, 2π]. There is something interesting going on
near the saddle. In Fig. C.2 we consider a narrower range of trajectory parameters,
and also zoom in closer to the saddle. A pair of caustics: envelopes of intersections of
nearby trajectories and the cusp at which they meet is clearly seen. Trajectory pattern
for yet a smaller range of trajectory parameters is shown in Fig. C.3. Incidentally, an
obvious but dramatic fact that trajectory parameters that differ very little can lead
to drastically different end points is illustrated in Fig. C.4; I believe that the actions
will also differ significantly. As stated above, the shadow region can be reached by
first exiting the basin, but the action of a trajectory that reaches the separatrix on the
shadow side, should generally be larger then the action of a trajectory on the other
side, which reaches the separatrix directly, without having to leave the basin first. Let
us make another zoom in the range of trajectory parameters: this is shown in Fig. C.5.
We see that the projection unto the (Q, P ) space of the most repelling eigenvector
in fact forms the cut, or the boundary of the shadow region. It must be noted that
if the plot was made over a larger range, the fan of trajectories actually intersects
the projection of the eigenvector. This may need a more thorough investigation: to
look at narrower and narrower ranges of trajectory parameters and observing if they
intersect the projection of the eigenvector closer and closer to the saddle. In any
case, it appears that the eigenvector defines the boundary of the shadow region only
infinitesimally close to the saddle.
There is another technicality. One of the checks that is to repeat calculations
at ever decreasing step sizes until the results do not change significantly (the other
check of course is the comparison with available analytics). We made the calculations
shown in Fig. C.5 at 1/10 the step size and aside from a slight shift of the trajectory
parameter of the most likely escape trajectory, nothing much changes qualitatively.

C.2

Caustics emanating from saddles

There are situations when something interesting appears as we continue to zoom in
to narrower and narrower range of trajectory parameters around the one that hits

147
the saddle. This is illustrated in Fig. C.6, where the parameters were η = 0.09 and
F ≈ 0.343 is here). A new type of caustic is apparent and does not seem to be a
direct continuation of the cusp pair formed by the trajectories in the larger parameter
range. For some other values of parameters, it was observed that as the parameter
range is varied from the relatively wide one to a narrower one, the position at which
the trajectories appear to intersect moves in closer to the saddle. Evidently there
is a caustic emanating from the saddle. In fact, caustics emanating from saddles
have been predicted in [31]. That paper argues that whenever the ratio of the larger
repelling eigenvalue to the smaller repelling eigenvalue λ = λlg /λsm is less then 3/2,
there is a caustic that emanates from the saddle point tangentially to the cut. In the
above example, λ = 1.37, our observation is consistent with this prediction. So far we
have not seen violations of this prediction. For example, when η = 0.13 and F ≈ 0.19
the behavior just described does not happen, and indeed, λ ≈ 1.92 there. However, a
more systematic numerical study is needed to verify this claim of [31], and we intend
to do this in the near future.

C.3

Rotation of eigenvectors

Shadows appear to exist also in the case of up → down transitions. In some cases this
phenomenon has been responsible for the seeming change in the slope of the effective
barrier height versus the driving strength. Upon the examination of trajectories it
was noticed that before the change in slope (say for F1 ) the most likely exit trajectory
approaches the saddle from one side and for after the change in slope (say for F2 > F1 )
it approaches the saddle from the other side. This lead me to suppose that there is
something akin to the first order transition, where there are actually several (locally
minimizing) paths that lead to the saddle, and the values of their actions switch.
Another possibility was simply a rotation of the eigenvector that defines the cut.
This turned out to be the case, and moreover, the calculations show that in general
this rotation is smooth, thus eliminating other exotic possibilities.
Fig. C.7 shows the snapshots of the fan of trajectories for three different values of

148
F for up → down transitions. All of these cases are actually very far from the up →
down bifurcation point. The probable reason that the slope of the effective barrier
height versus the driving strength appears to change is that it becomes harder to
hit the saddle, as is evident from the middle plot in Fig. C.7. Indeed, upon a very
painstaking bracketing, the slope change disappears, and nothing special is seen in
the function of the effective barrier height versus the driving strength.

149

1.5

0.5

−0.5

−1

−0.5

0.5

Figure C.1: A set of 50 trajectories that escape the attracting fixed point are shown
in blue, the most likely escape trajectory (which we know must pass through the
saddle) is in red and the relaxational trajectory from the saddle to the attracting
fixed point appears in black. Note, that these two are not just related by reversing
the sign of the friction, as would be in the case of detailed balance. The separatrix
of the lower-amplitude basin is shown in green.

150

−0.45

−0.5

−0.55

−0.6

−0.65

−0.7

−0.75

−0.8
0.45

0.5

0.55

0.6

0.65

0.7

0.75

0.8

Figure C.2: Same features as in Fig. C.1, but for a narrower range of parameters,
[4.5, 5.985] and closer to the saddle. The cusp and the pair of caustics are clearly
seen.

151

−0.7

−0.71

−0.72

−0.73

−0.74

−0.75

−0.76

−0.77

−0.78

−0.79

−0.8
0.5

0.52

0.54

0.56

0.58

0.6

0.62

Figure C.3: A set of 50 trajectories in the range of parameters [5.6328, 5.63398]. The
the region shown is even smaller then in Fig. C.2.

152

0.2

−0.2

−0.4

−0.6

−0.8

−0.4

−0.2

0.2

0.4

0.6

0.8

Figure C.4: Same features as in Fig. C.3, but showing the larger region of the (Q, P )
space. Some trajectories hit the separatrix very close to the saddle and some, for
a slightly different trajectory parameter, lead to a completely different region of the
(Q, P ) space.

153

−0.787
−0.7875
−0.788
−0.7885
−0.789

−0.7895
−0.79
−0.7905
−0.791
−0.7915
−0.792
0.508 0.5085 0.509 0.5095 0.51 0.5105 0.511 0.5115 0.512 0.5125 0.513

Figure C.5: A set of 50 trajectories in a very narrow window of trajectory parameters, [5.633489536, 5.63348954776]. The projection of the eigen-direction of the mostrepelling eigenvector is also shown in magenta.

154

0.7885

0.7890

0.7895

-0.2135

-0.2140

-0.2145

Figure C.6: Two sets of trajectories, one for a wide parameter region and one for a
much narrower one. The saddle point is not shown.

155

Figure C.7: Pattern of escape trajectories for up → down transitions at η = 0.15.
Top to bottom: patterns of escape trajectories at F corresponding to 0.8%,0.9% and
0.989998% of the hysteresis respectively. The separatrix (green) and the saddle point
(dark blue dot) are also shown.

156

Appendix D
Full Locking Criterion
Looking for fixed points of a set of m oscillators is an interesting problem by itself,
so what follows is a brief summary of this problem.
The dynamical equations for phase differences, φi = θi − θi−1 are:
φ̇−1 = ∆ω−1 + K0 cos φ0 − 2K−1 sin φ−1 + K−2 sin φ−2
φ̇0 = ∆ω0 + K1 sin φ1 + K−1 sin φ−1

φ̇1 = ∆ω1 + K2 sin φ2 − 2K1 sin φ1 − K0 cos φ0
φ̇2 = ∆ω2 + K3 sin φ3 − 2K2 sin φ2 + K1 sin φ1
...
φ̇m−1 = ∆ωm−1 + Km sin φm − 2Km−1 sin φm−1 + Km−2 sin φm−2
φ̇m = ∆ωm + Km+1 cos φm+1 − 2Km sin φm + Km−1 sin φm−1

φ̇m+1 = ∆ωm+1 + Km+2 sin φm+2 + Km sin φm
φ̇m+2 = ∆ωm+2 + Km+3 sin φm+3 − 2Km+2 sin φm+2 + Km+1 cos φm+1
Now isolate oscillators θ0 through θm from the rest of the chain, by setting K0 =
Km+1 = 0. The EOM are:

157

φ̇1

 
∆ω1

 
 
 
 
 φ̇2   ∆ω2  
 
 
 
 
 φ̇3  =  ∆ω3 +
 
 
 
 
 ...   ...  
 
 
φ̇m
∆ωm


−2K1

...

−2K2

...

−2K3 K4

...

Km−1

...

sin φ1



...   sin φ2 


...   sin φ3 


...
... 

−2Km
sin φm

The FP is:
−

−1 
−2K1

...

−2K2

...

−2K3 K4

...

sin φ1

∆ω1

 
 
 ∆ω2   sin φ2 
 
 
 ∆ω3  =  sin φ3 
 
 
 ...   ... 
 
sin φm
∆ωm

Km−1 −2Km

So there exist 2m FPs if all the entries in the vector on the lhs are ≤ 1. In a work by
Ermentrout and Kopell [94], it was proven that that if this is true, then one of the
2m possible FP is a sink. Therefore, as long as the entries of the lhs are ≤ 1, there
exists a global attractor.
We proceed and find the inverse of the square matrix. With the help of Mathematica I found the following patten:

−1 
m+1
 2
 Km−1
Km

m−1
K1
2(m−1)
K2

m−2
K1
2(m−2)
K2
3(m−2)
K3

m−3
K1
2(m−3)
K2
3(m−3)
K3

...

2·3
K2

2·2
K2

...

3·3
K3

3·2
K3

...


∆ω1

sin φ1

 

 

  ∆ω2   sin φ2 
 

 
3 
sin φ3 
∆ω3
K3  
=
 

...
...
...  
 
 

...   ∆ωm−1   sin φm−1 
sin φm
∆ωm
Km
(D.1)
K2

the inverted matrix is symmetric about both diagonals.
We now derive the probability distribution for the value of each row. The following

158
point of view is taken. Let us imagine an ensemble of frequency layouts for a given
layout of couplings. This will allows us to derive Pj (s, Kj ) - the probability of a value
in the jth row for a given Kj . Then, we imagine an ensemble of couplings and obtain
Pj (s) = Pj (s, K)G(K) dK. Here G(K) denotes the distribution of Ks (distribution
of ∆s is denoted by g(∆)).

Each row has the form Sj =

i=1 cj,i ∆i where Kj

distribution Pj (s, Kj ) is in principle given by

is included in these cs. The

g(∆1 , ∆2 , ...)δ(Sj − s) d∆1 d∆2 ..., but

this is a dead-end definition. However, its easy to calculate the characteristic function
of this Pj first, and then inverse Fourier-Transform it. Define the characteristic funcR
tion Φj = eikSj (∆1 ,∆2 ,...)g(∆1 , ∆2 , ...) d∆1 d∆2 .... Its easy to see that this is indeed

the characteristic function because upon expanding the exponential, we would obtain
the series in which each coefficient contains the nth moment. Now, because all the
∆s are IID, this Φj becomes Φj = eikcj,1 ∆1 g(∆1 ) d∆1 eikcj,2 ∆2 g(∆2 ) d∆2 .... So we
need to calculate

Z ∞

ikc∆

g(∆) d∆ =

−∞

Z ∞

eikc∆

−∞

λ/π
λ2 + ∆2

We can see (from the contour integral for example) that the answer is e−kcλ. Hence,
the characteristic function for jth row is
Φj = e−kλΣj ,
where Σj =

−kλ
, then e−kλΣj
i=1 cj,i So, since the characteristic function of g(∆) is e

is a characteristic function of
Pj (s, Kj ) =

(λΣj ) /π
(λΣj )2 + s2

We see from the matrix above that Σj is given by
[j (1 + 2 + ... + (m − j)) + (m − j + 1) (1 + 2 + ... + j)] = j(m−j+1)
. We should
Kj (m+1)
2Kj

now compute Pj (s) ≡

Pj (s, K)G(K) dK. This is not easy. However, what we are

after is the probability that a cluster of size m does not break. This means that |s| is
R1
less then 1 in every bond. In other words, this probability is given by Πm
j=1 −1 Pj (s) ds.

159
Let us take the s-integral first, and then deal with the K-integral. So, exchanging the
order of integration this way, the probability that a cluster of size m does not break
is
µZ ∞

¶ µZ ∞
Z 1
(λΣm ) /π
(λΣ1 ) /π
dKG(K)
ds ...
dKG(K)
ds
−1 (λΣ1 ) + s
−1 (λΣm ) + s
Z 1

2K
Z ∞ arctan
λj(m−j+1)
4µ
dK 
j=1
µ ¶m Y
m µZ ∞
³ µ´
arctan (aj x)
dx = f m,
π2
1 + x2
j=1

¶m Y

where aj = µλ j(m−j+1)
. This is a bit complicated, so let us first bound it: obtain lower

and upper limits. Note that the coefficient aj is smallest in the middle of the cluster
is smaller for smaller a
and biggest at the edges. Moreover, the area under arctan(ax)
1+x2
(and approaches ≈ 2.47 as a → ∞). So the lower bound can be obtained by setting
j = m/2 and the upper bound by setting all j = 1 for example. So, the upper limit
is:

u(m, µ/λ) =

π2

¶m ÃZ ∞

¢ !m
arctan µλ m2 x
dx
1 + x2

and the lower limit is
´ m
Z ∞ arctan µ 2 x
λ m(m+2)
l(m, µ/λ) =
dx
π2
1 + x2
¢ !m
µ ¶m ÃZ ∞
arctan µλ m22 x
dx
π2
1 + x2

¶m

where in the last term it was assumed that m À 1. Both expressions are functions
of µ/λ and m. We plot these functions below.

160

uHΜΛL, lHΜΛL
1.0
0.8
0.6
0.4
0.2

0.1

100

1000

ΜΛ

Figure D.1: The functions u(m, µ/λ) and l(m, µ/λ) plotted versus µ/λ for three values
of m: 1 (solid line), 3 (short dashed curve) and 7 (long dashed curve). In each case,
the upper and lower curve corresponds, respectively to functions u and l. Notice the
logarithmic scale of the x-axis.

161

Bibliography
[1] R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford: Oxford University
Press, 2001.
[2] E. M. Lifshitz and L. P. Pitayevskii, Physical Kinetics, Oxford: Pergamon Press,
1981.
[3] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
[4] H.W.C. Postma, I. Kozinsky, A. Husain, and M. L. Roukes, Appl. Phys. Lett.
86, 223105 (2005).
[5] R. Lifshitz and M. C. Cross, Nonlinear Dynamics of Nanomechanical and Micromechanical Resonators in Reviews of Nonlinear Dynamics and Complexity 1,
H. G. Schuster, editor, Weinhein: Wiley, 2008.
[6] I. Kozinsky, H. W. C. Postma, O. Kogan, A. Husain, and M. L. Roukes. Phys.
Rev. Lett. 99, 207201 (2007).
[7] M. Roukes, Phys. World 14, 25 (2001).
[8] M. LaHaye, O. Buu, B. Camarota, and K. Schwab, Science 304, 74 (2004).
[9] I. Katz, A. Retzker, R. Straub, and R. Lifshitz, Phys. Rev. Lett. 99, 040404
(2007).
[10] I. Katz, R. Lifshitz, A. Retzker, and R. Straub, New J. Phys. 10, 125023 (2008).
[11] D. Wahyu-Utami, and A. A. Clerk Phys. Rev. A, 78, 042323 (2008).

162
[12] I. Kozinsky, Nonlinear Nanoelectromechanical Systems, Ph.D. Thesis, Caltech,
(2007).
[13] H. A. Kramers, Physica VII, 4, 284 (1940).
[14] J. Kurkijarvi, Phys. Rev. B 6, 832 (1972).
[15] R. H. Victora, Phys. Rev. Lett. 63, 457 (1989).
[16] M. I. Dykman, and M. A. Krivoglaz, Physica 104A, 480-494 (1980).
[17] O. A. Tretiakov, Phys. Rev. B 67, 073303 (2003).
[18] K. A. Wiesenfeld, and E. Knobloch, Phys. Rev. A 26, 2946 (1982).
[19] R. S. Maier and D. L. Stein, Phys. Rev. Lett. 86, 3942 (2001).
[20] D. Ryvkine, M. I. Dykman, and B. Golding, Phys. Rev. E 69, 061102 (2004).
[21] M. I. Dykman, I. B. Schwartz, and M. Shapiro, Phys. Rev. E 72 021102 (2005).
[22] J. S. Aldrige and A. N. Cleland, Phys. Rev. Lett. 94, 156403 (2005).
[23] R. Almog, S. Zaitsev, O. Shtempluck, and E. Buks, Appl. Phys. Lett. 90, 013508
(2007).
[24] C. A. Stambaugh and H. B. Chan, Phys. Rev. B 73, 172302 (2006).
[25] A.N. Cleland and M.L. Roukes, J. Appl. Phys. 92, 2758 (2002).
[26] I. Siddiqi, R. Vijay, F. Pierre, C. M. Wilson, M. Metcalfe, C. Rigetti, L. Frunzio,
and M. H. Devoret, Phys. Rev. Lett. 93, 207002 (2004).
[27] Kihwan Kim, Myoung-Sun Heo, Ki-Hwan Lee, Hyoun-Jee Ha, Kiyoub Jang,
Heung-Ryoul Noh, and Wonho Jhe, Phys. Rev. A 72, 053402 (2005).
[28] B. Abdo, E. Segev, O. Shtempluck, and E. Buks, J. App. Phys. 101, 083909
(2007).

163
[29] R. Graham and H. Haken, Z. Physik 243, 289-302 (1971).
[30] R. Graham and T. Tél, Phys. Rev. A 31, 1109 (1985).
[31] M. I. Dykman, M. M. Millonas, V. N. Smelyanskiy, Phys. Lett. A 195, 53-58
(1994).
[32] R. Graham and T. Tél, Phys. Rev. A 33, 1322 (1986).
[33] A. Pikovsky, Synchronization: A Universal Concept in Nonlinear Science, New
York: Cambridge University Press, 2001.
[34] F. Varela, J. P. Lachaux, E. Rodriguez, J. Martinerie, Nature Rev. Neurosci. 2,
229 (2001).
[35] A. K. Engel and W. Singer, Trends Cogn. Sci. 5, 16 (2001).
[36] M. Rabinovich, P. Varona, A. Selverston, and H. Abarbanel, Rev. Mod. Phys.
78, 1213 (2006).
[37] K. Wiesenfeld, Physica B 222, 315 (1996).
[38] B. R. Trees, V. Saranathan, and D. Stroud, Phys. Rev. E 71, 016215 (2005).
[39] B. R. Trees, S. Natu, and D. Stroud, Phys. Rev. B 72, 214524 (2005).
[40] A. G. Vladimirov, G. Kozyreff, and P. Mandel, Europhys. Lett. 61, 613 (2003).
[41] J. F. Heagy, T. L. Carroll, and L. M. Pecora, Phys. Rev. E 50, 1874 (1994).
[42] S. K. Esser, S. L. Hill, and G. Tononi, Sleep 30, 1617 (2007).
[43] J. J. Collins, and I. N. Stewart, J. Nonlinear Sci. 3, 349 (1993).
[44] L. Glass, Nature 410, 277 (2001).
[45] J. A. Acerbón, L. L. Bonilla, C. J. Pérez-Vicente, F. Ritort, and R. Spigler, Rev.
Mod. Phys. 77, 137-185 (2005).

164
[46] Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, 39, 420, H. Araki, editor, New York:
Springer, 1975.
[47] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Berlin: Springer,
1984.
[48] B. C. Daniels, S. T. M. Dissanayake, and B. R. Trees, Phys. Rev. E 67, 026216
(2003).
[49] F. Rogister, and R. Roy, Phys. Rev. Lett. 98, 104101 (2007).
[50] N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group,
Reading: Addison-Wesley, 1992.
[51] S. H. Strogatz, and R. E. Mirollo, J. Phys. A: Math. Gen. 21, L699-L705 (1988).
[52] S. H. Strogatz, R. E. Mirollo, Physica D 31, 143-168 (1988).
[53] M. I. Dykman, and M. A. Krivoglaz, Theory of Nonlinear Oscillator Interacting
with a Medium in Soviet Physics Reviews, 5, 265-441, I. M. Khalatnikov, editor,
New York: Harwood Academic, 1984.
[54] M. I. Freidlin, and A. D. Wentzell, Random perturbations of dynamical systems,
New York: Springer, 1998.
[55] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry, and
the Natural Sciences, Berlin: Springer-Verlag, 2004.
[56] P. J. Holmes, and D. A. Rand, Journal of Sound and Vibration, 44, 237-253
(1976).
[57] P. Holmes, and D. Rand, International Journal of Non-linear Mechanics, 15,
449-458 (1980).
[58] A. H. Nayfeh, and D. T. Mook, Nonlinear oscillations, New York: John Wiley,
1979.

165
[59] A. H. Nayfeh, Perturbation methods, New York: John Wiley & Sons, 2000.
[60] N. N. Bogoliubov, and Y. A. Mitropolsky, Asymptotic Methods in the Theory of
Nonlinear Oscillations, New York: Gordon and Breach Science Publishers, 1961.
[61] F. Verhulst, Nonlinear Differential Equations and Dynamical Systems, Berlin:
Springer-Verlag, 1990.
[62] L. D. Landau, and E. M. Lifshitz, Mechanics, Oxford: Pergamon Press, 1989.
[63] O. Kogan, Phys. Rev. E 76, 037203 (2007).
[64] R. K. Pathria, Statistical Mechanics, Oxford: Butterworth-Heinemann, 1996.
[65] M.I. Dykman, and M. A. Krivoglaz, Sov. Phys. JETP 50, 30 (1979).
[66] R. Graham, Z. Physik B 26 281-290 (1977).
[67] M. I. Dykman, Phys. Rev. A. 42, 2020 (1990).
[68] M. I. Dykman in Stochastic and Chaotic Dynamics in the Lakes, D. S. Broomhead, E. A. Luchinskaya, P. V. E. McClintock, and T. Mullin, editors, New York:
Melville, 1999.
[69] M. I. Dykman, and V. N. Smelyanskiy, Superlattices and Microstructures 23,
495 (1998).
[70] R. Phythian, J. Phys. A: Math. Gen. 10, 777 (1977).
[71] B. Jouvert, and R. Phythian, Phys. Rev. A 19, 1350 (1979).
[72] R. P. Feynman, and A. R. Hibbs, Quantum Mehcanics and Path Integrals, New
York: McGraw-Hill, 1965.
[73] L. S. Schulman, Techniques and Applications of Path Integration, Mineola:
Dover, 2005.
[74] D. L. Stein, and R. S. Maier, SIAM J. Appl. Math. 57, 752 (1997).

166
[75] V. P. Maslov, and M. V. Fedoriuk, Semiclassical Approximations in Quanutm
Mechanics, Dordrecht: Reidel, 1981.
[76] V. I. Arnold, Mathematical Methods of Classical Mechanics, Berlin: Springer,
1978.
[77] A.P. Dmitriev, M.I. D’yakonov, Zh. Eksp. Teor. Fiz. 90, 1430-1439, (1986); also
in Sov. Phys. JETP 63, 838 (1986).
[78] J. Guckenheimer, and P. Holmes, Nonlinear Oscillators, Dynamical Systems and
Bifurcations of Vector Fields, New York: Springer-Verlag, 1987.
[79] V. Arnold, Ordinary Differential Equations, 3rd edition, New York: Springer,
1992.
[80] O. Kogan, arXiv:0805.0972v2. A shorter version co-authored with M. I. Dykman
is to be submitted to Phys. Rev. E.
[81] V. A. Chinarov, M. I. Dykman, and V. N. Smelyanskiy, Phys. Rev. E 47, 2448
(1993).
[82] R. S. Maier, and D. L. Stein, Phys. Rev. Lett. 85, 1358 (2000).
[83] H. Whitney, Ann. Math. 62, 374 (1955).
[84] V. I. Arnold, Catastrophe theory, Berlin: Springer, 1984.
[85] H. R. Jauslin, Physica A 144, 179 (1987).
[86] M. V. Day, Stochasitcs 20, 121 (1987).
[87] R. Graham and T. Tél, Phys. Rev. Lett. 52, 9 (1984).
[88] V. N. Smelyanskiy, M. I. Dykman, and R. S. Maier, Phys. Rev. E 55, 2369
(1996).
[89] O. Kogan, and M. I. Dykman, to be submitted.

167
[90] M. Bennett, M. F. Schatz, H. Rockwood, and K. Wiesenfeld, Proc. Roy. Soc.
Series A 458, 563 (2002).
[91] A. T. Winfree, J. Theor. Biol. 16, 15 (1967).
[92] M. C. Cross, J. L. Rogers, R. Lifshitz, and A. Zumdieck, Phys. Rev. E 73, 036205
(2006).
[93] S. H. Strogatz, Physica D 143, 1-20 (2000).
[94] G. Ermentrout, and N. Kopell, Siam J. Math. Anal. 15, 215 (1984).
[95] S. K. Ma, C. Dasgupta, and C. K. Hu, Phys. Rev. Lett. 43, 1434 (1979).
[96] C. Dasgupta, and S. K. Ma, Phys. Rev. B 22, 1305 (1980).
[97] D. S. Fisher, Phys. Rev. B 50, 3799 (1994).
[98] D. S. Fisher, and A. P. Young, Phys. Rev. B 58, 9131 (1998).
[99] E. Altman, Y. Kafri, A. Polkovnikov, and G. Refael, Phys. Rev. Lett. 93, 150402
(2004).
[100] E. Altman, Y. Kafri, A. Polkovnikov, and G. Refael, Phys. Rev. Lett. 100,
170402 (2008).
[101] O. Kogan, J. L. Rogers, M. C. Cross, G. Refael, submitted to Phys. Rev. E.;
also on arXiv:0810.3075v2.
[102] K. Huang, Statistical Mechanics, New York: Wiley, 1963.
[103] T. Lee, G. Refael, M. C. Cross, O. Kogan, and J. L. Rogers, to appear.
[104] H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, Prog. Theor. Phys. 77, 1005
(1987).
[105] H. Daido, Phy. Rev. Lett. 61, 231 (1988).